Topfolio
Tutorial

Python Exploratory Data Analysis (EDA): The Complete Step-by-Step Workflow

Master the complete 6-stage Python Exploratory Data Analysis (EDA) framework. Learn structured data inspection, missing value imputation, IQR outlier detection, distribution analysis, correlation heatmaps, and feature profiling.

Anuj SainiAug 24, 202613 min read

When junior analysts receive a messy dataset, they often jump straight into training machine learning models or building visualizations without inspecting data types, nullity mechanisms, or distributional skew. The result? Biased metrics, silent pipeline failures, and flawed business conclusions.

Senior analysts follow a structured, repeatable EDA playbook. Whether you are analyzing e-commerce transactions, user engagement metrics, or medical clinical trials, this 6-stage workflow guarantees that you uncover data quality bugs, understand feature relationships, and extract actionable insights.

For foundational Python syntax before diving into this playbook, review our Python Pandas Data Analysis Guide or follow the comprehensive Data Analyst Career Track.

Interactive Practice Environment

Sharpen your data cleaning, filtering, and statistical analysis skills with instant browser-based execution on Topfolio Practice.

Stage 1: Structural Health & Data Type Audit

Before running calculations, you must understand your dataset's dimensionality, column dtypes, memory utilization, and summary statistics.

1.1 Environment Setup & Configuration

Configure Pandas and Seaborn to display all columns and produce clean, legible chart aesthetics:

python
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import warnings
 
# Display configuration
pd.set_option('display.max_columns', None)
pd.set_option('display.max_rows', 100)
pd.set_option('display.float_format', lambda x: f'{x:.3f}')
sns.set_theme(style='whitegrid', palette='muted')
warnings.filterwarnings('ignore')
 
print("EDA Environment Initialized.")

1.2 Ingesting and Inspecting Structural Metadata

Let's generate a realistic, messy dataset with numerical signals, categorical groups, skewed features, and injected nulls:

python
from sklearn.datasets import make_classification
 
# Generate synthetic tabular dataset
np.random.seed(42)
X, y = make_classification(
    n_samples=1200,
    n_features=8,
    n_informative=5,
    n_redundant=2,
    n_classes=2,
    random_state=42
)
 
feature_names = [f'feature_{i}' for i in range(8)]
df = pd.DataFrame(X, columns=feature_names)
df['target'] = y
 
# Inject real-world categorical columns and messiness
df['region'] = np.random.choice(['North', 'South', 'East', 'West'], size=len(df), p=[0.4, 0.3, 0.2, 0.1])
df['customer_tier'] = np.random.choice(['Bronze', 'Silver', 'Gold', 'Platinum'], size=len(df), p=[0.5, 0.3, 0.15, 0.05])
 
# Inject missing values
df.loc[df.sample(frac=0.08, random_state=1).index, 'feature_0'] = np.nan
df.loc[df.sample(frac=0.05, random_state=2).index, 'feature_3'] = np.nan
df.loc[df.sample(frac=0.04, random_state=3).index, 'region'] = np.nan
 
# Inject a heavy right-skewed revenue feature
df['annual_spend'] = np.random.exponential(scale=5000, size=len(df)) + 500
 
print(f"Dataset Dimensions: {df.shape[0]} rows × {df.shape[1]} columns")

1.3 Structural Audit Commands

Run these three standard diagnostics immediately upon opening any dataset:

python
# 1. Inspect data types and non-null counts
df.info()
 
# 2. Summary statistics for continuous numerical columns
numerical_summary = df.describe().T
numerical_summary['skewness'] = df.select_dtypes(include=np.number).skew()
print("\nNumerical Feature Summary:")
print(numerical_summary[['count', 'mean', 'std', 'min', '50%', 'max', 'skewness']])
 
# 3. Summary statistics for categorical columns
categorical_summary = df.describe(include=['object']).T
print("\nCategorical Feature Summary:")
print(categorical_summary)
ColumnNon-Null CountDtypeInterpretation
feature_0 to feature_71,104–1,200float64Continuous numerical features
target1,200int64Binary classification target (0/1)
region1,152objectNominal categorical feature (4 categories)
customer_tier1,200objectOrdinal categorical feature (4 levels)
annual_spend1,200float64Skewed financial metric

Stage 2: Missing Value Diagnostics & Imputation

Missing values are not simply "blank cells." In statistics, missing data falls into three distinct classifications:

  1. MCAR (Missing Completely at Random): Missingness is purely random and unrelated to any variable (e.g., random network packet drop).
  2. MAR (Missing at Random): Missingness depends on observed data, not the missing value itself (e.g., women being less likely to report weight in a survey, but recorded via gender column).
  3. MNAR (Missing Not at Random): Missingness directly relates to the unobserved value (e.g., high-income earners refusing to disclose salary).

2.1 Quantifying Missingness

python
def missing_data_report(data: pd.DataFrame) -> pd.DataFrame:
    missing_count = data.isnull().sum()
    missing_pct = (missing_count / len(data)) * 100
    report = pd.DataFrame({
        'Missing_Count': missing_count,
        'Missing_Percentage': missing_pct,
        'Dtype': data.dtypes
    })
    return report[report['Missing_Count'] > 0].sort_values(by='Missing_Count', ascending=False)
 
print(missing_data_report(df))

2.2 Visualizing Missingness Patterns

A nullity heatmap displays whether missing values occur simultaneously across multiple columns:

python
plt.figure(figsize=(10, 5))
sns.heatmap(df.isnull(), cbar=False, cmap='magma', yticklabels=False)
plt.title("Missing Value Matrix (White bars indicate missingness)", fontsize=13, pad=12)
plt.xlabel("Columns")
plt.tight_layout()
plt.show()

2.3 Strategic Imputation Implementation

Never drop missing rows indiscriminately. Apply tailored imputation strategies based on data type and distribution:

python
df_cleaned = df.copy()
 
# Strategy 1: Impute skewed numerical features with Median
# (Mean is vulnerable to extreme values; median preserves central tendency)
numerical_cols_with_nulls = ['feature_0', 'feature_3']
for col in numerical_cols_with_nulls:
    median_val = df_cleaned[col].median()
    # Add a binary missingness indicator column to preserve signal for downstream ML
    df_cleaned[f'{col}_was_missing'] = df_cleaned[col].isnull().astype(int)
    df_cleaned[col].fillna(median_val, inplace=True)
 
# Strategy 2: Impute categorical features with Mode or 'Unknown' category
df_cleaned['region'].fillna(df_cleaned['region'].mode()[0], inplace=True)
 
# Verify zero remaining nulls
print(f"Remaining Missing Values: {df_cleaned.isnull().sum().sum()}")

Stage 3: Outlier Detection & Anomaly Treatment

Outliers can represent measurement errors, fraudulent transactions, or rare organic events. In data analytics, two primary methods isolate numerical anomalies:

  1. IQR Rule (Interquartile Range): Non-parametric, robust against non-normal distributions.
  2. Z-Score Method: Parametric, best suited for Gaussian (normal) distributions.
                    IQR Method Outlier Boundaries
       Lower Bound                                Upper Bound
    Q1 - 1.5 * IQR     Q1      Median      Q3    Q3 + 1.5 * IQR
 [--- Outliers ---]----|=========|=========|----[--- Outliers ---]
                       |<------- IQR ----->|

3.1 Implementing the IQR Outlier Detector

python
def detect_outliers_iqr(data: pd.DataFrame, column: str) -> tuple[pd.DataFrame, float, float]:
    q25 = data[column].quantile(0.25)
    q75 = data[column].quantile(0.75)
    iqr = q75 - q25
    
    lower_bound = q25 - (1.5 * iqr)
    upper_bound = q75 + (1.5 * iqr)
    
    outliers = data[(data[column] < lower_bound) | (data[column] > upper_bound)]
    return outliers, lower_bound, upper_bound
 
outliers, lower_b, upper_b = detect_outliers_iqr(df_cleaned, 'annual_spend')
print(f"Outlier boundaries for 'annual_spend': [{lower_b:.2f}, {upper_b:.2f}]")
print(f"Detected {len(outliers)} outliers ({len(outliers) / len(df_cleaned) * 100:.2f}% of total rows)")

3.2 Z-Score Outlier Detection

python
from scipy import stats
 
z_scores = np.abs(stats.zscore(df_cleaned['feature_1']))
z_outliers = df_cleaned[z_scores > 3.0]
print(f"Features with |Z-score| > 3 in feature_1: {len(z_outliers)}")

3.3 Outlier Remediation: Winsorization (Clipping)

Instead of dropping rows (which deletes valuable data in other columns), winsorize (clip) extreme outliers to the 1st and 99th percentiles:

python
def winsorize_column(series: pd.Series, lower_percentile=0.01, upper_percentile=0.99) -> pd.Series:
    lower_cap = series.quantile(lower_percentile)
    upper_cap = series.quantile(upper_percentile)
    return series.clip(lower=lower_cap, upper=upper_cap)
 
# Clip extreme annual spend without dropping customer transactions
df_cleaned['annual_spend_clipped'] = winsorize_column(df_cleaned['annual_spend'])
 
print(f"Original Max Spend: ${df_cleaned['annual_spend'].max():.2f}")
print(f"Clipped Max Spend:  ${df_cleaned['annual_spend_clipped'].max():.2f}")

Stage 4: Univariate Distribution Analysis

Univariate analysis inspects features individually to evaluate variance, skewness, spread, and modal peaks.

4.1 Numerical Features: Histograms + KDE

python
features_to_plot = ['feature_0', 'feature_1', 'feature_2', 'annual_spend']
 
fig, axes = plt.subplots(2, 2, figsize=(14, 9))
axes = axes.flatten()
 
for idx, col in enumerate(features_to_plot):
    sns.histplot(df_cleaned[col], kde=True, ax=axes[idx], color='teal', bins=30)
    skew = df_cleaned[col].skew()
    axes[idx].set_title(f"{col} (Skew: {skew:.2f})", fontsize=12, fontweight='bold')
    axes[idx].set_xlabel("Value")
    axes[idx].set_ylabel("Frequency")
 
plt.tight_layout()
plt.show()

4.2 Transforming Skewed Distributions

Highly skewed variables (like annual_spend with skew > 1.5) destabilize statistical regressions. Apply a Log Transformation ($y = \ln(x + 1)$):

python
df_cleaned['log_annual_spend'] = np.log1p(df_cleaned['annual_spend_clipped'])
 
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))
sns.histplot(df_cleaned['annual_spend_clipped'], kde=True, ax=ax1, color='salmon')
ax1.set_title(f"Before Log: Raw Spend (Skew: {df_cleaned['annual_spend_clipped'].skew():.2f})")
 
sns.histplot(df_cleaned['log_annual_spend'], kde=True, ax=ax2, color='mediumseagreen')
ax2.set_title(f"After Log: Transformed Spend (Skew: {df_cleaned['log_annual_spend'].skew():.2f})")
plt.tight_layout()
plt.show()

4.3 Categorical Frequencies

python
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 4.5))
 
order_tier = ['Bronze', 'Silver', 'Gold', 'Platinum']
sns.countplot(data=df_cleaned, x='customer_tier', order=order_tier, ax=ax1, palette='Blues_r')
ax1.set_title("Customer Tier Distribution", fontweight='bold')
ax1.set_xlabel("Tier")
 
sns.countplot(data=df_cleaned, x='region', ax=ax2, palette='Set2')
ax2.set_title("Regional Distribution", fontweight='bold')
ax2.set_xlabel("Region")
 
plt.tight_layout()
plt.show()

Stage 5: Bivariate & Interaction Analysis

Bivariate analysis explores dependencies between features and the target outcome.

5.1 Numerical Feature vs. Binary Target (Boxplots)

python
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
 
features = ['feature_0', 'feature_1', 'log_annual_spend']
for idx, feature in enumerate(features):
    sns.boxplot(x='target', y=feature, data=df_cleaned, ax=axes[idx], palette='coolwarm')
    axes[idx].set_title(f"{feature} by Target Class", fontweight='bold')
    axes[idx].set_xlabel("Target (0 = Churned, 1 = Retained)")
 
plt.tight_layout()
plt.show()

5.2 Categorical vs. Target (Normalized Crosstabs)

To evaluate if certain categories convert at higher rates, generate a 100% stacked bar chart using normalized frequency tables:

python
# Compute proportion of target class per region
region_crosstab = pd.crosstab(df_cleaned['region'], df_cleaned['target'], normalize='index') * 100
 
print("Conversion / Retention Rate by Region (%):")
print(region_crosstab)
 
ax = region_crosstab.plot(kind='bar', stacked=True, figsize=(9, 5), color=['#e74c3c', '#2ecc71'])
plt.title("Target Proportion by Geographic Region", fontsize=13, fontweight='bold')
plt.xlabel("Region")
plt.ylabel("Percentage (%)")
plt.legend(['Class 0', 'Class 1'], title='Target', bbox_to_anchor=(1.02, 1), loc='upper left')
plt.xticks(rotation=0)
plt.tight_layout()
plt.show()

Stage 6: Correlation Profiling & Feature Engineering

Correlation heatmaps reveal collinearity (features that duplicate information) and highlight the strongest predictors of the target.

6.1 Masked Correlation Heatmap

python
# Compute correlation matrix for numeric columns
numeric_df = df_cleaned.select_dtypes(include=np.number).drop(columns=['target'])
corr_matrix = numeric_df.corr(method='pearson')
 
# Generate an upper triangle mask to prevent visual redundancy
mask = np.triu(np.ones_like(corr_matrix, dtype=bool))
 
plt.figure(figsize=(12, 8))
sns.heatmap(
    corr_matrix,
    mask=mask,
    annot=True,
    fmt='.2f',
    cmap='vlag',
    vmin=-1.0,
    vmax=1.0,
    linewidths=0.75,
    cbar_kws={'label': 'Pearson Correlation Coefficient'}
)
plt.title("Feature Collinearity Matrix (Lower Triangle)", fontsize=14, fontweight='bold', pad=14)
plt.tight_layout()
plt.show()
Correlation Interpretation Matrix:
| Value Range | Correlation Strength | Action Required |
| :--- | :--- | :--- |
| 0.80 to 1.00 | Very Strong Collinearity | Drop one feature or apply PCA to prevent variance inflation |
| 0.40 to 0.79 | Moderate Correlation | Excellent candidate for interaction feature engineering |
| -0.10 to 0.10 | Negligible Linear Relation | Test non-linear relationships or binning |

6.2 Feature Engineering: Interactions, Binning & One-Hot Encoding

Translate EDA insights into machine-learning-ready features:

python
# 1. Interaction Feature (Ratio between two correlated features)
df_cleaned['feature_0_1_ratio'] = df_cleaned['feature_0'] / (np.abs(df_cleaned['feature_1']) + 1e-5)
 
# 2. Quantile Binning (Transforming continuous spend into 4 balanced tiers)
df_cleaned['spend_quartile'] = pd.qcut(
    df_cleaned['annual_spend_clipped'],
    q=4,
    labels=['Low_Spend', 'Mid_Low_Spend', 'Mid_High_Spend', 'High_Spend']
)
 
# 3. One-Hot Encoding nominal and binned categorical columns
df_encoded = pd.get_dummies(
    df_cleaned,
    columns=['region', 'customer_tier', 'spend_quartile'],
    drop_first=True
)
 
print(f"Final Encoded Dataset Shape: {df_encoded.shape[0]} rows × {df_encoded.shape[1]} columns")
df_encoded.head(3)

Complete Reproducible EDA Script

You can copy and execute this entire workflow in any Jupyter Notebook or script:

python
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
from scipy import stats
 
def run_automated_eda(filepath: str, target_col: str = None):
    """
    Executes a structured 6-stage EDA diagnostic on any CSV dataset.
    """
    print("=" * 60)
    print("TOPFOLIO AUTOMATED EDA DIAGNOSTIC ENGINE")
    print("=" * 60)
    
    # 1. Load Data
    data = pd.read_csv(filepath)
    print(f"[STAGE 1] Ingested {data.shape[0]} rows and {data.shape[1]} columns.")
    
    # 2. Missing Values Check
    nulls = data.isnull().sum()
    null_cols = nulls[nulls > 0]
    if not null_cols.empty:
        print("\n[STAGE 2] Missing Columns Identified:")
        for col, count in null_cols.items():
            print(f"  - {col}: {count} nulls ({count/len(data)*100:.1f}%)")
    else:
        print("\n[STAGE 2] No missing values detected.")
        
    # 3. Outlier Audit via IQR
    print("\n[STAGE 3] Numerical Outlier Audit (IQR Method):")
    num_cols = data.select_dtypes(include=np.number).columns
    for col in num_cols:
        q25, q75 = data[col].quantile(0.25), data[col].quantile(0.75)
        iqr = q75 - q25
        outliers = data[(data[col] < q25 - 1.5*iqr) | (data[col] > q75 + 1.5*iqr)]
        if len(outliers) > 0:
            print(f"  - {col}: {len(outliers)} outliers ({len(outliers)/len(data)*100:.1f}%)")
            
    # 4. Target Relationship (if specified)
    if target_col and target_col in data.columns:
        print(f"\n[STAGE 4 & 5] Target Correlations with '{target_col}':")
        corrs = data[num_cols].corr()[target_col].sort_values(ascending=False)
        print(corrs)
        
    print("\n" + "=" * 60)
    print("EDA DIAGNOSTIC COMPLETE")
    print("=" * 60)
    return data

Summary & What to Learn Next

Exploratory Data Analysis is not a passive inspection step — it is an active risk-management and hypothesis-generating discipline. By systematically auditing data types, managing missingness mechanisms, handling extreme values via winsorization, and visualizing collinear relationships, you ensure that every downstream model or executive dashboard rests on pristine data.

Next Steps in Your Analytics Roadmap:

Frequently Asked Questions

What is the difference between Exploratory Data Analysis (EDA) and Feature Engineering?

EDA is the discovery phase where you understand variable distributions, identify anomalies, inspect missingness mechanisms, and validate assumptions. Feature engineering is the transformation phase where you create new predictive indicators, encode categorical features, scale numerical ranges, or bin values based on insights gained during EDA.

When should I use Median imputation instead of Mean imputation?

Use Median imputation for skewed numerical distributions or variables with extreme outliers (such as salaries, transaction amounts, or order sizes), because the mean is heavily distorted by extreme values. Use Mean imputation only when numerical data follows a symmetric, normal distribution.

How does the Interquartile Range (IQR) rule detect outliers in Pandas?

The IQR rule computes the difference between the 75th percentile (Q3) and 25th percentile (Q1): IQR = Q3 - Q1. Any data point falling below Q1 - 1.5 * IQR or above Q3 + 1.5 * IQR is flagged as a statistical outlier. This method is non-parametric and does not assume a normal distribution.

What is the difference between Pearson and Spearman correlation in heatmaps?

Pearson correlation measures linear relationships between continuous variables and assumes normality. Spearman correlation measures monotonic relationships (whether variables increase together regardless of linearity) using rank orders, making it robust against non-linear patterns and extreme outliers.

How can I prevent data leakage during EDA and data cleaning?

Never calculate global statistics (like dataset-wide means, medians, or scaling parameters) across your full dataset before splitting. Always split into training and testing sets first, fit your imputers and scalers strictly on the training set, and apply those fitted values to the test set.

Anuj Saini

Written by

Anuj SainiFounder & Lead Instructor

Founder at Topfolio with 6+ years in data & analytics across JPMC, Ultrahuman, and high-growth startups. Sat on hiring panels, reviewed 500+ resumes, and writes practical SQL & data guides.

View all articles

Related Articles

Tutorial

Customer Analytics in Python: Cohort Analysis, RFM Segmentation & LTV

Master customer analytics in Python with Pandas and Seaborn. Calculate RFM scores, build customer segments, run monthly cohort retention heatmaps, and estimate Customer Lifetime Value (LTV).

13 min readRead
Tutorial

Market Basket Analysis in Python: Support, Confidence, Lift & Apriori Explained

Master Market Basket Analysis and Association Rule Mining in Python. Learn the mathematical intuition behind Support, Confidence, and Lift, one-hot encode transaction baskets, and implement Apriori and FP-Growth using Pandas and Mlxtend.

12 min readRead
Tutorial

REST APIs for Data Analysts in Python: Authentication, Pagination & JSON Normalization

Master REST API data extraction in Python. Learn how to handle Bearer tokens and API keys, loop through offset and cursor pagination, flatten nested JSON with pd.json_normalize(), and build fault-tolerant pipelines with automatic retries.

12 min readRead