Chapter 30: Exploratory Data Analysis

Getting Oriented: Data Inventory and Summary Statistics

The first ten minutes with any dataset should be the same every time: load it, count things, look at types, and ask a handful of blunt questions. How many rows and columns? What is each column's type — continuous numeric, integer count, nominal category, ordinal, timestamp, or free text? What fraction of values are missing in each column? What does the very first row look like, and the very last?

In pandas, this orientation pass costs almost nothing:

import pandas as pd
import numpy as np

df = pd.read_csv("titanic.csv")

print(df.shape)          # (891, 12)
print(df.dtypes)         # column types
print(df.head())
print(df.isnull().mean().sort_values(ascending=False))

The isnull().mean() call immediately reveals that cabin is missing for 77 % of passengers — a fact that shapes every analysis involving cabin assignment.

Descriptive Statistics for Continuous Variables

For numeric columns, df.describe() returns count, mean, standard deviation, and the five-number summary (min, Q1, median, Q3, max). This is enormously useful but easy to misread. Consider the standard deviation: it is sensitive to outliers in a way the interquartile range (IQR = Q3 − Q1) is not. When mean and median diverge substantially, suspect skewness.

Formal skewness and excess kurtosis help quantify what the eye sees in a histogram:

A right-skewed distribution (long tail to the right) has positive skewness. Excess kurtosis above zero indicates heavier tails than a Gaussian, which matters for outlier sensitivity of mean-based methods.

from scipy import stats

numeric_cols = df.select_dtypes(include="number").columns
for col in numeric_cols:
    s = df[col].dropna()
    print(f"{col:20s}  mean={s.mean():.2f}  median={s.median():.2f}"
          f"  skew={stats.skew(s):.2f}  kurt={stats.kurtosis(s):.2f}")

For the Titanic fare column you will see a mean of roughly $32 but a median near $14, with skewness above 4 — a classic right-skewed count-like variable that benefits from log-transformation before modeling.

Descriptive Statistics for Categorical Variables

For object or categorical columns, df[col].value_counts() and df[col].nunique() are the go-to functions. Always check the cardinality (number of distinct values) early: a column with 850 unique values in a 891-row dataset is probably a near-identifier, not a useful feature.

for col in df.select_dtypes(include=["object", "category"]).columns:
    vc = df[col].value_counts(normalize=True)
    print(f"\n{col}  ({df[col].nunique()} unique)")
    print(vc.head(5).to_string())

Common Pitfalls

• Reading describe() uncritically. The count row shows how many non-null values exist, not the total rows. A column with a dramatically lower count than the others is missing data silently.
• Ignoring string columns. df.describe() by default omits object columns unless you pass include='all'. Missing that means missing an entire class of data quality issues.
• Treating ID columns as features. Columns with monotonically increasing integers or unique strings are not meaningful predictors; they inflate cardinality and can leak row order into models.

Univariate Analysis: Understanding Each Variable in Isolation

Before comparing any two variables, spend time with each one individually. Univariate analysis answers: What values does this column take? How are they distributed? Are there anomalies? The goal is to build a mental model of each feature so that downstream bivariate and multivariate findings can be interpreted correctly.

Continuous Variables: Histograms and KDE Plots

A histogram divides the range of a continuous variable into bins and counts observations per bin. The choice of bin width is critical — too few bins and you miss structure; too many and you see noise. A simple rule is Sturges' rule ($k = \lceil \log_2 n \rceil + 1$), but Freedman-Diaconis ($h = 2 \cdot \text{IQR}(x) \cdot n^{-1/3}$) is more robust to outliers. In practice, overlaying a kernel density estimate (KDE) on the histogram combines the histogram's bin-count accuracy with a smooth density curve:

import matplotlib.pyplot as plt
import seaborn as sns

fig, axes = plt.subplots(1, 2, figsize=(12, 4))

# Histogram + KDE for fare
sns.histplot(df["Fare"].dropna(), kde=True, bins=40, ax=axes[0])
axes[0].set_title("Raw Fare Distribution")
axes[0].set_xlabel("Fare (£)")

# Log-transformed fare looks much more Gaussian
sns.histplot(np.log1p(df["Fare"].dropna()), kde=True, bins=40, ax=axes[1])
axes[1].set_title("Log(1 + Fare) Distribution")
axes[1].set_xlabel("log(1 + Fare)")

plt.tight_layout()
plt.savefig("fare_distribution.png", dpi=150)

The log-transformed version is closer to symmetric — confirming the earlier skewness finding and suggesting log-fare is the better representation for linear models.

Box Plots: Five-Number Summary at a Glance

A box plot encodes the five-number summary visually: the box spans Q1 to Q3 (the IQR), the central line is the median, and the whiskers extend to $Q1 - 1.5 \cdot \text{IQR}$ and $Q3 + 1.5 \cdot \text{IQR}$. Points beyond the whiskers are plotted individually and conventionally flagged as potential outliers — though this is a heuristic, not a verdict.

fig, ax = plt.subplots(figsize=(6, 4))
sns.boxplot(y=df["Age"].dropna(), ax=ax, color="steelblue")
ax.set_title("Age Distribution")
plt.tight_layout()

For age, the box plot reveals a median around 28, mild right skew (whisker longer on top), and a handful of points near 70-80 that the IQR rule flags — but which are plausible ages, not errors.

Categorical Variables: Bar Charts and Frequency Tables

For categorical variables, a bar chart of normalized counts conveys the same information as value<em>counts(normalize=True), but visually:

fig, axes = plt.subplots(1, 2, figsize=(12, 4))

# Passenger class
df["Pclass"].value_counts().sort_index().plot.bar(ax=axes[0], rot=0, color="steelblue")
axes[0].set_title("Passenger Class Counts")
axes[0].set_xlabel("Class")

# Embarked
df["Embarked"].value_counts().plot.bar(ax=axes[1], rot=0, color="coral")
axes[1].set_title("Port of Embarkation")

plt.tight_layout()

Third class is the most common (55 %), and Southampton is overwhelmingly the dominant embarkation port (72 %). Both are facts you need before interpreting any class-stratified or port-stratified finding.

When to Transform

The decision to transform a variable during EDA is about understanding, not cosmetics:

• Log or log1p for right-skewed positive quantities (income, price, count variables). Compresses the upper tail; makes multiplicative relationships additive.
• Square root for count data with moderate skew (softer than log).
• Rank or quantile normalization when the shape is irregular and you want rank-order preserved without distributional assumptions.
• Standardization (z-score) to compare variables on different scales — but remember it does not change shape, only location and scale.

Always plot both raw and transformed distributions side by side to confirm the transformation achieves its purpose.

Checking for Modes and Multimodality

A unimodal distribution has one peak; a bimodal distribution has two. Bimodality often signals a mixture of subpopulations. In the Titanic age column, a subtle secondary peak near age 1-5 reflects infant passengers — which might warrant separate treatment or a derived is</em>child feature. KDE plots with bandwidth tuning are the best tool for detecting modes without binning artifacts:

fig, ax = plt.subplots(figsize=(8, 4))
for bw in [0.5, 1.0, 2.0]:
    df["Age"].dropna().plot.kde(bw_method=bw, ax=ax, label=f"bw={bw}")
ax.legend()
ax.set_title("Age KDE with Different Bandwidths")

Larger bandwidth smooths away the infant peak; smaller bandwidth preserves it but introduces noise. The right bandwidth sits where structure stabilizes across a range of values.

Outlier Detection: Finding and Handling Anomalous Values

Outliers are data points that lie far from the bulk of the distribution. They matter for multiple reasons: they can represent genuine extreme events (a whale customer, a record-breaking measurement), data entry errors (age = 999, price = -1), or the natural heavy-tail behavior of a phenomenon. EDA cannot decide which kind an outlier is — that requires domain knowledge — but it can find them systematically and present them for scrutiny.

The IQR Rule

The most widely used rule flags a value as a potential outlier if it lies more than 1.5 × IQR below Q1 or above Q3. This was proposed by Tukey (1977) and has the property that for a Gaussian distribution it catches roughly 0.7 % of observations — rare enough to be worth investigating.

def iqr_outliers(series: pd.Series) -> pd.Series:
    """Return boolean mask — True where IQR outlier."""
    q1, q3 = series.quantile([0.25, 0.75])
    iqr = q3 - q1
    return (series < q1 - 1.5 * iqr) | (series > q3 + 1.5 * iqr)

for col in ["Age", "Fare", "SibSp", "Parch"]:
    mask = iqr_outliers(df[col].dropna())
    pct = 100 * mask.mean()
    print(f"{col:10s}: {mask.sum():4d} outliers ({pct:.1f}%)")

For Fare, you will see roughly 17 % flagged — evidence that this rule's threshold is too aggressive for a highly skewed variable. Always interpret the percentage in light of the distribution shape; the IQR rule is designed for roughly symmetric data.

Z-Score Method

The z-score measures how many standard deviations a value is from the mean:

A common threshold is $|z| > 3$, catching about 0.3 % of a Gaussian distribution.

from scipy.stats import zscore

df["fare_z"] = zscore(df["Fare"].fillna(df["Fare"].median()))
outlier_rows = df[df["fare_z"].abs() > 3][["Name", "Fare", "Pclass", "fare_z"]]
print(outlier_rows.sort_values("fare_z", ascending=False).head(10))

The extreme-fare passengers are predominantly first-class travelers, with fares above £500 — this is a real phenomenon (expensive suites), not an error.

Modified Z-Score (Robust to Skewness)

Both the IQR rule and z-score can be misleading for skewed distributions because the mean and standard deviation are themselves pulled by the outliers. A more robust alternative uses the median absolute deviation (MAD):

Values with $|M_i| > 3.5$ are flagged. The constant 0.6745 makes MAD consistent with the standard deviation under normality.

def modified_zscore(series: pd.Series) -> pd.Series:
    med = series.median()
    mad = (series - med).abs().median()
    return 0.6745 * (series - med) / (mad + 1e-10)  # avoid div-by-zero

df["fare_mz"] = modified_zscore(df["Fare"].dropna())
print((df["fare_mz"].abs() > 3.5).sum(), "outliers via modified z-score")

Visualization: Box Plots and Strip Plots

Numerical flags are only the start — visualization reveals context. A strip plot (jittered scatter) overlaid on a box plot shows every point, distinguishing the near-outliers from the extreme ones:

fig, ax = plt.subplots(figsize=(8, 5))
sns.boxplot(x="Pclass", y="Fare", data=df, ax=ax,
            flierprops=dict(marker="o", color="red", markersize=4))
sns.stripplot(x="Pclass", y="Fare", data=df, ax=ax,
              alpha=0.3, jitter=True, color="steelblue", size=3)
ax.set_title("Fare by Passenger Class with Outliers Highlighted")
plt.tight_layout()

This immediately shows that high fares are almost exclusively in first class, and that the outliers cluster there rather than being scattered randomly — pointing to suite pricing rather than data error.

What to Do With Outliers

Once identified, outliers require a decision:

• Investigate first. Look up the rows. Does the value make sense in context?
• Domain clamp. If you know the valid range (e.g., age must be 0–120), values outside that range are errors to be corrected or nullified.
• Winsorize. Replace values above the 99th percentile with the 99th percentile value. This preserves the row while limiting influence. Use scipy.stats.mstats.winsorize.
• Transform. Log-transforming a right-skewed variable often makes former outliers look ordinary by compressing the upper tail.
• Separate stratum. Sometimes outliers form their own legitimate segment (e.g., enterprise customers) that deserves its own model.
• Never silently delete without documenting what was removed and why.

Bivariate Analysis: Relationships Between Two Variables

Once you understand each variable in isolation, the investigative question shifts: how do two variables relate? The answer depends on the types of the two variables.

Continuous vs. Continuous: Scatter Plots and Correlation

A scatter plot is the primary tool for two continuous variables. You are looking for:

• Direction: does one variable increase as the other increases (positive), decrease (negative), or neither?
• Strength: how tightly do points cluster around the trend line?
• Form: is the relationship linear, curved, or non-existent?
• Outliers: are there points far from the cloud that may be driving any apparent relationship?

fig, ax = plt.subplots(figsize=(7, 5))
sns.scatterplot(x="Age", y="Fare", hue="Survived", data=df,
                alpha=0.6, palette={0: "red", 1: "steelblue"}, ax=ax)
sns.regplot(x="Age", y="Fare", data=df, scatter=False,
            color="gray", line_kws={"linestyle": "--"}, ax=ax)
ax.set_title("Age vs. Fare, colored by Survival")

The regression line (from regplot) shows any linear trend; coloring by survival lets you check whether the relationship differs between groups — the beginning of multivariate thinking.

Pearson's correlation coefficient $r$ quantifies linear association on a $[-1, 1]$ scale:

$r = 0$ does not mean no relationship; a U-shaped relationship can have $r \approx 0$. Spearman's rank correlation $\rho$ is more appropriate when the relationship is monotone but not linear, or when outliers are present:

from scipy.stats import pearsonr, spearmanr

clean = df[["Age", "Fare"]].dropna()
r_p, p_p = pearsonr(clean["Age"], clean["Fare"])
r_s, p_s = spearmanr(clean["Age"], clean["Fare"])
print(f"Pearson  r={r_p:.3f}  p={p_p:.4f}")
print(f"Spearman r={r_s:.3f}  p={p_s:.4f}")

Correlation Matrices and Heatmaps

When you have more than two numeric variables, computing all pairwise correlations at once reveals which pairs deserve deeper investigation:

import seaborn as sns

corr = df[["Survived", "Pclass", "Age", "SibSp", "Parch", "Fare"]].corr(method="spearman")

fig, ax = plt.subplots(figsize=(7, 6))
sns.heatmap(corr, annot=True, fmt=".2f", cmap="coolwarm",
            center=0, square=True, ax=ax)
ax.set_title("Spearman Correlation Matrix")
plt.tight_layout()

The heatmap immediately surfaces: Survived correlates negatively with Pclass (higher class = lower number = higher survival) and positively with Fare — and those two correlations are themselves large because fare and class are collinear.

Continuous vs. Categorical: Box Plots and Violin Plots

When one variable is categorical and the other continuous, the question is: does the distribution of the continuous variable differ across groups?

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Box plot: median and spread by survival
sns.boxplot(x="Survived", y="Fare", data=df, palette="Set2", ax=axes[0])
axes[0].set_xticklabels(["Died", "Survived"])
axes[0].set_title("Fare by Survival Status")

# Violin plot: full distribution shape
sns.violinplot(x="Pclass", y="Age", data=df, palette="Set3",
               inner="quartile", ax=axes[1])
axes[1].set_title("Age Distribution by Passenger Class")

plt.tight_layout()

The violin plot reveals that third-class passengers tend to be slightly younger than first-class passengers, with a more irregular age distribution — including that infant peak.

Categorical vs. Categorical: Crosstabs and Grouped Bar Charts

For two categorical variables, a crosstab (contingency table) reveals co-occurrence counts, and a normalized version reveals conditional probabilities:

# Survival rate by sex
crosstab = pd.crosstab(df["Sex"], df["Survived"], normalize="index")
print(crosstab.round(3))

# Visual
crosstab.plot.bar(rot=0, figsize=(6, 4),
                  color=["salmon", "steelblue"],
                  title="Survival Rate by Sex")
plt.ylabel("Proportion")
plt.legend(["Died", "Survived"])
plt.tight_layout()

The output shows female survival at ~74 % versus male survival at ~19 % — one of the strongest single-variable signals in the dataset.

Testing for Association

Beyond visualizing, you can formally test whether an association is statistically significant. For categorical-categorical relationships, the chi-squared test of independence applies:

from scipy.stats import chi2_contingency

contingency = pd.crosstab(df["Sex"], df["Survived"])
chi2, p, dof, expected = chi2_contingency(contingency)
print(f"chi2={chi2:.2f}, p={p:.2e}, dof={dof}")

A p-value near zero confirms what the chart shows: sex and survival are not independent. But remember, statistical significance does not equal causal importance — it means the association is unlikely to be sampling noise.

Groupby, Aggregation, and Pivot Tables

Some of the most actionable insights in EDA come not from individual observations but from groups of them. Groupby-aggregation answers questions like: What is the average fare for each class? Does survival rate differ by embarkation port? How does age distribution change across decade-of-birth cohorts? These segmented summaries are the analytical backbone of business intelligence and the foundation of feature engineering.

The Split-Apply-Combine Pattern

The pandas groupby API implements the split-apply-combine paradigm:

1. Split the DataFrame into subgroups based on one or more keys.
2. Apply an aggregation, transformation, or filter to each group.
3. Combine the results back into a new DataFrame.

# Single-key aggregation
survival_by_class = (
    df.groupby("Pclass")["Survived"]
    .agg(["mean", "count", "sum"])
    .rename(columns={"mean": "survival_rate", "count": "n", "sum": "n_survived"})
    .round(3)
)
print(survival_by_class)

This confirms the expected pattern: first-class survival is 63 %, second-class 47 %, third-class 24 %.

Multiple Aggregations with Named Aggregation

Pandas' named aggregation syntax (agg(col=(column, func))) allows combining different statistics for different columns cleanly:

group_stats = (
    df.groupby(["Pclass", "Sex"])
    .agg(
        n                = ("PassengerId", "count"),
        survival_rate    = ("Survived",    "mean"),
        median_age       = ("Age",         "median"),
        median_fare      = ("Fare",        "median"),
        fare_p90         = ("Fare",        lambda x: x.quantile(0.9)),
    )
    .reset_index()
    .round(2)
)
print(group_stats.to_string(index=False))

The result is a compact segment profile: third-class males have a survival rate of 14 %, a median age of 25, and a median fare of £8. First-class females have a survival rate of 97 %.

Pivot Tables

A pivot table is a two-dimensional crosstab that summarizes a numeric value as a function of two categorical keys — the row key and the column key. In pandas:

pivot = df.pivot_table(
    values="Survived",
    index="Pclass",
    columns="Sex",
    aggfunc="mean",
    margins=True,
    margins_name="Overall",
).round(3)
print(pivot)

The output is a 4×3 table (three classes plus overall row, two sexes plus overall column) of survival rates. The margins=True argument appends row and column totals — here, marginal survival rates.

Pivot tables work equally well for non-binary outcomes:

fare_pivot = df.pivot_table(
    values="Fare",
    index="Pclass",
    columns="Embarked",
    aggfunc="median",
    margins=True,
).round(1)
print(fare_pivot)

Custom Aggregation Functions

When built-in aggregation functions are insufficient, pass a callable:

def coef_of_variation(x):
    """Standard deviation divided by mean — relative variability."""
    return x.std() / x.mean() if x.mean() != 0 else np.nan

df.groupby("Pclass")["Fare"].agg(coef_of_variation).rename("CV_fare")

Coefficient of variation reveals that fare variability is highest in first class (CV ≈ 0.85) — evidence that first class had a range of cabin types with very different prices.

Visualizing Group Comparisons

Aggregate summaries become most persuasive when plotted:

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Survival rate by class and sex (grouped bar)
group_stats.pivot(index="Pclass", columns="Sex", values="survival_rate").plot.bar(
    ax=axes[0], rot=0, color=["steelblue", "salmon"]
)
axes[0].set_title("Survival Rate by Class and Sex")
axes[0].set_ylabel("Rate")
axes[0].legend(title="Sex")

# Median fare by class and sex (grouped bar)
group_stats.pivot(index="Pclass", columns="Sex", values="median_fare").plot.bar(
    ax=axes[1], rot=0, color=["steelblue", "salmon"]
)
axes[1].set_title("Median Fare by Class and Sex")
axes[1].set_ylabel("Fare (£)")

plt.tight_layout()

Rolling and Expanding Aggregations for Time Series

Groupby is not limited to categorical splits. When your data has a time dimension, rolling windows reveal trends:

# Hypothetical: daily sales DataFrame
# sales_df["revenue"].rolling(window=7, min_periods=1).mean()  # 7-day moving average
# sales_df["revenue"].expanding().mean()                       # cumulative mean

These techniques apply the same split-apply-combine logic over time windows rather than discrete categories.

Pitfalls

• Simpson's paradox. A trend observed in grouped data can reverse when the groups are combined. Always check both disaggregated and aggregate views. The Titanic is a textbook example: overall, passengers in third class had lower fares AND lower survival, but within each sex the class gradient is consistent.
• Unequal group sizes. A mean of 100 % survival from a group of 2 is less meaningful than 90 % from a group of 100. Always include group counts in your aggregation output.
• Multiple groupby keys and memory. Grouping on high-cardinality keys produces many small groups; aggregation results can explode in size. Check df.groupby(keys).ngroups before aggregating on unknown cardinality.

Multivariate Analysis and Hypothesis-Driven EDA

Univariate and bivariate analyses reveal individual variables and pairwise relationships. But real datasets contain interactions: the effect of variable A on outcome Y may depend on the level of variable B. Multivariate EDA explores these higher-order structures while maintaining the discipline of asking explicit questions and checking them against evidence.

The Hypothesis-Driven Mindset

Effective EDA is not aimless chart generation — it is a cycle of hypothesize → visualize → quantify → refine. Before each plot, write down the question it answers. After each plot, decide: does the evidence support the hypothesis, refute it, or suggest a new one? This discipline keeps your analysis focused and produces a narrative that stakeholders can follow.

For the Titanic, a natural initial hypothesis is:

H1: Passenger class determines survival, primarily through access to lifeboats.

We already know marginal survival rates by class. The next question is whether class explains the sex effect, or whether both sex and class have independent effects.

Faceted Plots: Conditioning on a Third Variable

Faceted plots (also called small multiples) condition the bivariate relationship on a third categorical variable, making interactions visible. In seaborn, FacetGrid or the col/row parameters of most plot functions handle this:

g = sns.FacetGrid(df.dropna(subset=["Age"]),
                  col="Pclass", row="Survived",
                  height=3, aspect=1.2,
                  margin_titles=True)
g.map(sns.histplot, "Age", bins=20, kde=True, color="steelblue")
g.set_axis_labels("Age", "Count")
g.set_titles(row_template="Survived={row_name}",
             col_template="Class {col_name}")
g.figure.suptitle("Age Distribution by Class and Survival", y=1.02)
plt.savefig("faceted_age.png", dpi=150, bbox_inches="tight")

With 6 panels (3 classes × 2 survival outcomes), you can now see: among third-class passengers who died, the distribution peaks in the 20s. Among those who survived, the infant peak is more prominent — consistent with a "children first" policy applied unevenly across classes.

ggplot2: Facets and Layers in R

ggplot2's grammar of graphics makes multivariate layering especially expressive. The facet<em>grid and facet</em>wrap functions apply the same faceting logic:

library(tidyverse)

titanic <- read_csv("https://raw.githubusercontent.com/datasciencedojo/datasets/master/titanic.csv") |>
  mutate(
    Survived = factor(Survived, labels = c("Died", "Survived")),
    Pclass   = factor(Pclass, labels = c("1st","2nd","3rd"))
  )

ggplot(titanic |> drop_na(Age), aes(x = Age, fill = Sex)) +
  geom_histogram(aes(y = after_stat(density)), bins = 20,
                 alpha = 0.6, position = "identity") +
  geom_density(aes(color = Sex), linewidth = 0.8, fill = NA) +
  facet_grid(Survived ~ Pclass) +
  scale_fill_manual(values  = c("female"="steelblue", "male"="salmon")) +
  scale_color_manual(values = c("female"="steelblue", "male"="salmon")) +
  labs(title = "Age by Sex, Class, and Survival Outcome",
       x = "Age", y = "Density") +
  theme_minimal(base_size = 11)

This six-panel plot encodes four variables simultaneously: age (x), density (y), sex (color), class (column), survival (row).

Pair Plots with a Grouping Variable

A pair plot (scatter matrix) shows all pairwise scatter plots and univariate distributions together. The diagonal typically shows KDE or histogram; the off-diagonal shows scatter plots:

cols = ["Age", "Fare", "Pclass", "SibSp", "Survived"]
subset = df[cols].dropna().copy()
subset["Survived"] = subset["Survived"].map({0: "Died", 1: "Survived"})

g = sns.pairplot(subset, hue="Survived",
                 palette={"Died": "salmon", "Survived": "steelblue"},
                 diag_kind="kde", corner=True,
                 plot_kws={"alpha": 0.35, "s": 15})
g.figure.suptitle("Pair Plot by Survival", y=1.02)

Correlation Does Not Imply Causation — and Neither Does Any EDA

EDA is generative, not confirmatory. Every finding is a hypothesis for further investigation, not a conclusion. The strong association between Sex == female and Survived == 1 does not by itself prove that being female caused survival — it is consistent with a lifeboat-loading policy, but also with wealth differences (women more concentrated in first class) or other confounders. Testing causality requires experimental design or causal inference methods that lie beyond EDA.

Building the EDA Narrative

By the end of a thorough EDA, you should be able to write a short memo with:

• Data quality findings: missing value patterns, potential errors, outliers requiring treatment.
• Distributional summaries: which variables are skewed, which have outliers, which need transformation.
• Key relationships: the strongest predictors of the outcome variable, notable interactions, surprising null findings.
• Hypotheses carried forward: questions that EDA raised but cannot answer — e.g., whether class and sex interact in a statistically significant way, which requires a formal model.

This memo becomes the justification for every preprocessing and feature engineering decision you make next.