Exploratory Data Analysis (EDA): A Step-by-Step Framework

Every failed model traces back to a data problem the analyst never saw. Exploratory data analysis (EDA) is the discipline of interrogating a dataset before you model it: its shape, its distributions, its gaps, its relationships, and its lies. This framework gives you a repeatable seven-step workflow, complete with copy-pasteable Python, so you can move from raw CSV to modeling-ready data without guessing.

Quick Takeaways

• EDA is a structured sequence: load, profile, clean, univariate, bivariate, multivariate, then feature decisions.
• Run df.info(), df.describe(), and df.isna().sum() before any chart. They catch most problems in seconds.
• Use skewness, IQR fences, and correlation matrices to drive preprocessing choices, not intuition.
• Document every decision (imputation, capping, dropping) so the pipeline is reproducible.

What Is Exploratory Data Analysis?

Exploratory data analysis is the process of summarizing a dataset’s main characteristics, typically with statistics and visualization, to uncover structure, anomalies, and hypotheses. John Tukey formalized the approach in the 1970s. His principle still holds: look at the data before you test anything.

EDA answers five practical questions:

1. What does each column represent, and is its type correct?
2. How much data is missing, duplicated, or impossible?
3. How is each variable distributed?
4. Which variables relate to the target and to each other?
5. What transformations does the data need before modeling?

Skipping EDA does not save time. It moves the debugging to after training, where it costs more.

Step 1: Load and Inspect the Dataset

Start with structure. Check row count, column count, dtypes, and memory footprint.

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

# Configure display so wide dataframes are not truncated
pd.set_option("display.max_columns", 50)
sns.set_theme(style="whitegrid")

# Load data; parse dates up front to avoid object-typed timestamps
df = pd.read_csv("customers.csv", parse_dates=["signup_date"])

print(df.shape)            # (rows, columns)
print(df.head())           # eyeball the first records
df.info(memory_usage="deep")  # dtypes, non-null counts, memory

Look for three red flags immediately:

• Numeric columns stored as object: usually currency symbols, commas, or “N/A” strings.
• Dates stored as object: breaks time-based features.
• ID columns stored as numeric: they should never be averaged or correlated.

Define the Target and Variable Types

Separate columns by role. This drives every later choice.

target = "churned"

# Auto-split by dtype; review manually afterward
num_cols = df.select_dtypes(include="number").columns.drop(target).tolist()
cat_cols = df.select_dtypes(include=["object", "category"]).columns.tolist()

print(f"Numeric: {len(num_cols)} | Categorical: {len(cat_cols)}")

Step 2: Profile Data Quality

Quality profiling quantifies how broken the data is. Measure missingness, duplicates, and cardinality.

# Missing values as count and percentage, sorted worst-first
missing = (
    df.isna().sum()
    .to_frame("missing_count")
    .assign(missing_pct=lambda x: (x["missing_count"] / len(df) * 100).round(2))
    .query("missing_count > 0")
    .sort_values("missing_pct", ascending=False)
)
print(missing)

# Exact duplicate rows
print("Duplicates:", df.duplicated().sum())

# Cardinality: high-cardinality categoricals need special encoding
print(df[cat_cols].nunique().sort_values(ascending=False))

Classify the Missingness Mechanism

How you impute depends on why values are missing.

Test for pattern by comparing groups:

# Does missingness in 'income' differ by churn status?
df["income_missing"] = df["income"].isna().astype(int)
print(df.groupby(target)["income_missing"].mean())

A large gap between groups signals MAR or MNAR. Keep the indicator as a feature.

Step 3: Clean and Standardize

Apply fixes only after profiling. Record each one.

# Remove exact duplicates
df = df.drop_duplicates()

# Fix a numeric column polluted by currency formatting
df["monthly_spend"] = (
    df["monthly_spend"]
    .astype(str)
    .str.replace(r"[$,]", "", regex=True)
    .replace("nan", np.nan)
    .astype(float)
)

# Standardize inconsistent category labels
df["plan"] = df["plan"].str.strip().str.lower()

# Impute: median for skewed numerics, mode for categoricals
df["monthly_spend"] = df["monthly_spend"].fillna(df["monthly_spend"].median())
df["plan"] = df["plan"].fillna(df["plan"].mode()[0])

Use the median over the mean when a distribution is skewed. The mean chases outliers; the median does not.

Step 4: Univariate Analysis

Study one variable at a time. The goal is shape, center, spread, and tails.

Numeric Variables

# Summary statistics plus skewness and kurtosis
summary = df[num_cols].describe().T
summary["skew"] = df[num_cols].skew()
summary["kurtosis"] = df[num_cols].kurt()
print(summary.round(2))

# Histogram + boxplot grid
fig, axes = plt.subplots(len(num_cols), 2, figsize=(12, 3 * len(num_cols)))
for i, col in enumerate(num_cols):
    sns.histplot(df[col], kde=True, ax=axes[i, 0])
    sns.boxplot(x=df[col], ax=axes[i, 1])
plt.tight_layout()
plt.show()

Interpret skewness with these thresholds:

• |skew| < 0.5: approximately symmetric.
• 0.5 to 1: moderate skew; consider a transform for linear models.
• > 1: heavy skew; apply log1p or a Box-Cox transform.

Detect Outliers

The IQR method flags values beyond 1.5 times the interquartile range from the quartiles.

def iqr_outliers(series: pd.Series, k: float = 1.5) -> pd.Series:
    """Return a boolean mask of IQR-based outliers."""
    q1, q3 = series.quantile([0.25, 0.75])
    iqr = q3 - q1
    return (series < q1 - k * iqr) | (series > q3 + k * iqr)

outlier_report = {c: iqr_outliers(df[c]).mean() * 100 for c in num_cols}
print(pd.Series(outlier_report).round(2).sort_values(ascending=False))

Never delete outliers by reflex. Verify whether each is an error (age of 240) or a legitimate extreme (a whale customer). Cap with clip() for errors; keep or flag legitimate values.

Categorical Variables

for col in cat_cols:
    print(df[col].value_counts(normalize=True).head(10).round(3))
    print("---")

# Visualize class balance of the target
sns.countplot(x=target, data=df)
plt.show()
print(df[target].value_counts(normalize=True))

Check the target balance now. A 95/5 split means accuracy is a misleading metric, and you will need AUC-PR, F1, or resampling.

Step 5: Bivariate Analysis

Bivariate analysis measures how each feature relates to the target and to other features. Match the method to the variable types.

from scipy import stats

# Numeric vs. binary target: compare distributions
sns.boxplot(x=target, y="monthly_spend", data=df)
plt.show()

# Mann-Whitney U: non-parametric test of group difference
g0 = df.loc[df[target] == 0, "monthly_spend"]
g1 = df.loc[df[target] == 1, "monthly_spend"]
stat, p = stats.mannwhitneyu(g0, g1)
print(f"Mann-Whitney U p-value: {p:.4g}")

# Categorical vs. target: churn rate per category
print(df.groupby("plan")[target].agg(["mean", "count"]))

# Chi-square test of independence
ct = pd.crosstab(df["plan"], df[target])
chi2, p, dof, _ = stats.chi2_contingency(ct)
print(f"Chi-square p-value: {p:.4g}")

Prefer Spearman correlation when relationships are monotonic but nonlinear, or when outliers remain. Pearson only captures linear association.

Step 6: Multivariate Analysis

Multivariate analysis reveals redundancy and interactions that pairwise views miss.

Correlation Heatmap

corr = df[num_cols + [target]].corr(method="spearman")

# Mask the upper triangle to halve the visual noise
mask = np.triu(np.ones_like(corr, dtype=bool))
plt.figure(figsize=(10, 8))
sns.heatmap(corr, mask=mask, annot=True, fmt=".2f", cmap="coolwarm", center=0)
plt.show()

Flag any feature pair with |correlation| above 0.85. Keep one, or combine them.

Multicollinearity with VIF

The Variance Inflation Factor (VIF) quantifies how well other predictors explain a feature. It matters most for linear and logistic regression.

from statsmodels.stats.outliers_influence import variance_inflation_factor
from statsmodels.tools.tools import add_constant

X = add_constant(df[num_cols].dropna())
vif = pd.DataFrame({
    "feature": X.columns,
    "VIF": [variance_inflation_factor(X.values, i) for i in range(X.shape[1])],
})
print(vif[vif["feature"] != "const"].sort_values("VIF", ascending=False))

A VIF above 5 warrants attention. Above 10, drop or merge the feature. Tree-based models tolerate collinearity far better, but it still dilutes feature importance.

Dimensionality Check with PCA

from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA

# PCA is scale-sensitive: always standardize first
X_scaled = StandardScaler().fit_transform(df[num_cols].dropna())
pca = PCA().fit(X_scaled)

# Cumulative explained variance: how many components capture 90%?
cum_var = np.cumsum(pca.explained_variance_ratio_)
print("Components for 90% variance:", np.argmax(cum_var >= 0.90) + 1)

If a few components explain most variance, your features carry heavy redundancy.

Step 7: Make Feature Engineering Decisions

EDA ends with decisions. Translate findings into a preprocessing plan.

Standardization vs. Normalization

Wrap decisions in a ColumnTransformer and Pipeline to prevent data leakage:

from sklearn.compose import ColumnTransformer
from sklearn.pipeline import Pipeline
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler, OneHotEncoder
from sklearn.model_selection import train_test_split

X = df.drop(columns=[target])
y = df[target]

# Stratify to preserve class balance in both splits
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, stratify=y, random_state=42
)

numeric_pipe = Pipeline([
    ("impute", SimpleImputer(strategy="median")),
    ("scale", StandardScaler()),
])
categorical_pipe = Pipeline([
    ("impute", SimpleImputer(strategy="most_frequent")),
    ("encode", OneHotEncoder(handle_unknown="ignore")),
])

preprocessor = ColumnTransformer([
    ("num", numeric_pipe, num_cols),
    ("cat", categorical_pipe, cat_cols),
])
# Fit on train only; transform both. This prevents leakage.

Split before you fit any transformer. Fitting a scaler on the full dataset leaks test-set statistics into training.

Automating EDA at Scale

Manual EDA does not scale to hundreds of columns. Use profiling libraries for a first pass, then investigate flagged columns by hand.

from ydata_profiling import ProfileReport

# minimal=True speeds up wide or large dataframes
report = ProfileReport(df, title="Churn EDA", minimal=True)
report.to_file("eda_report.html")

For datasets beyond memory, sample intelligently:

# Stratified sample keeps target distribution intact
sample = df.groupby(target, group_keys=False).apply(
    lambda x: x.sample(frac=0.1, random_state=42)
)

Real-World EDA Use Cases

Predicting Customer Churn

EDA on a telecom dataset typically exposes three signals: month-to-month contracts churn far more than annual ones, tenure correlates negatively with churn, and monthly charges are right-skewed. Decisions follow directly: one-hot encode contract type, log-transform charges, and use stratified splitting because churn is imbalanced.

Handling Missing Data in Real-Time Sensor Streams

Industrial sensors drop readings in bursts, not randomly. Plot missingness against time with msno.matrix() (missingno library). If gaps cluster during maintenance windows, the mechanism is MAR. Use forward-fill with a limit or interpolation within short gaps, and flag long gaps rather than imputing them.

Fraud Detection

Fraud rates sit below 1%. EDA reveals that transaction amount is heavy-tailed and that fraud concentrates at unusual hours. Prioritize Spearman correlations, RobustScaler, and time-derived features such as hour-of-day. Evaluate with AUC-PR, not accuracy.

Retail Demand Forecasting

Decompose the series to expose trend and seasonality before choosing a model:

from statsmodels.tsa.seasonal import seasonal_decompose

# Weekly seasonality in daily sales data
result = seasonal_decompose(sales["units"], model="additive", period=7)
result.plot()
plt.show()

Common EDA Mistakes to Avoid

• Exploring the test set. Do EDA on training data only, or findings leak into modeling.
• Trusting averages. Always pair the mean with spread and a distribution plot.
• Confusing correlation with causation. A strong r supports a hypothesis, not a conclusion.
• Ignoring Simpson’s paradox. A trend can reverse when you segment by a lurking variable. Always check subgroups.
• Dropping missing rows by default. Listwise deletion biases results under MAR and MNAR.
• Skipping documentation. Log each cleaning decision so another analyst can reproduce it.

Frequently Asked Questions

What are the main steps of exploratory data analysis?

The main steps are loading and inspecting the data, profiling quality (missing values, duplicates), cleaning and standardizing, univariate analysis, bivariate analysis, multivariate analysis, and making feature engineering decisions. Each step feeds the next and produces documented findings.

What is the difference between EDA and data cleaning?

Data cleaning fixes specific errors: wrong types, duplicates, and impossible values. EDA is broader. It includes cleaning, but also examines distributions and relationships to guide modeling choices. Cleaning is a step inside EDA.

Which Python libraries are best for EDA?

Pandas handles manipulation and summary statistics. Matplotlib and Seaborn handle visualization. SciPy and statsmodels provide statistical tests. ydata-profiling automates first-pass reports, and PySpark handles distributed data.

How do you handle outliers during EDA?

Detect them with IQR fences, modified Z-scores, or Isolation Forest. Then investigate their cause. Correct or cap errors, keep legitimate extremes, and consider robust methods like RobustScaler or tree-based models when outliers are real.