Exploratory Data Analysis

Good EDA is not about running every possible plot. It is about answering specific questions before you model. For this dataset the questions are: What drives house prices? Is the target clean? Are there outliers that will break linear models? Does geography matter?

Setup and First Look

from sklearn.datasets import fetch_california_housing
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

housing = fetch_california_housing(as_frame=True)
df = housing.frame

print(df.shape)      # (20640, 9)
print(df.head(3))
#    MedInc  HouseAge  AveRooms  AveBedrms  Population  AveOccup  Latitude  Longitude  MedHouseVal
# 0    8.33      41.0      6.98       1.02       322.0      2.56     37.88    -122.23         4.53
# 1    8.30      21.0      6.24       0.97      2401.0      2.11     37.86    -122.22         3.59
# 2    7.26      52.0      8.29       1.07       496.0      2.80     37.85    -122.24         3.46

print(df.isnull().sum())
# All zeros — no missing values

Tip

Always check .dtypes and .isnull().sum() before anything else. Knowing the data is all-numeric with no nulls tells you that you do not need imputation or encoding — you can spend your time on feature engineering instead.

Target Distribution

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

# Raw distribution
axes[0].hist(df['MedHouseVal'], bins=60, color='steelblue', edgecolor='none', alpha=0.8)
axes[0].axvline(5.0, color='red', linestyle='--', linewidth=1.5, label='Ceiling at 5.0')
axes[0].set_xlabel('MedHouseVal ($100k)')
axes[0].set_ylabel('Count')
axes[0].set_title('Target Distribution (Raw)')
axes[0].legend()

# Log-transformed distribution
log_target = np.log1p(df['MedHouseVal'])
axes[1].hist(log_target, bins=60, color='teal', edgecolor='none', alpha=0.8)
axes[1].set_xlabel('log(1 + MedHouseVal)')
axes[1].set_title('Target Distribution (Log-transformed)')

plt.tight_layout()
plt.show()

print(df['MedHouseVal'].skew().round(3))
# 1.079  — moderately right-skewed

print(np.log1p(df['MedHouseVal']).skew().round(3))
# 0.313  — much more symmetric after log transform

Warning

Notice the spike at 5.0. That vertical bar at the ceiling is not natural variation — it is the data collection cap. A normal distribution does not have a wall at one end. Any evaluation metric you compute includes predictions for those 965 capped blocks, which means your error metrics for the highest-value neighborhoods are artificially compressed.

Should You Log-Transform the Target?

Log-transforming the target (y = log1p(MedHouseVal)) often improves linear regression performance on skewed targets. The tradeoff is that your model now predicts in log space, and you must exponentiate predictions back before reporting to stakeholders.

# Log-transform and inverse-transform pattern
log_target = np.log1p(df['MedHouseVal'])

# After fitting a model and getting predictions in log space:
# log_predictions = model.predict(X_test)
# predictions_in_original_units = np.expm1(log_predictions)

Tree-based models (Random Forest, Gradient Boosting) do not need this transformation — they are not sensitive to target skew. Linear models benefit from it.

Correlation Analysis

corr = df.corr()

plt.figure(figsize=(9, 7))
mask = np.triu(np.ones_like(corr, dtype=bool))
sns.heatmap(
    corr,
    mask=mask,
    annot=True,
    fmt='.2f',
    cmap='coolwarm',
    center=0,
    square=True,
    linewidths=0.5
)
plt.title('Correlation Matrix — California Housing')
plt.tight_layout()
plt.show()

# Correlations with target, sorted
print(corr['MedHouseVal'].sort_values(ascending=False))
# MedHouseVal    1.000
# MedInc         0.688   ← dominant predictor
# AveRooms       0.151
# HouseAge       0.106
# AveBedrms     -0.047
# AveOccup      -0.023
# Population    -0.025
# Longitude     -0.046
# Latitude      -0.144

Success

MedInc has a correlation of 0.69 with the target. That is unusually strong for a single feature. Income predicts house prices better than house age, room count, or location alone. Your model should reflect this — feature importance from tree models and coefficients from linear models should both place MedInc at the top.

Income vs. House Value

plt.figure(figsize=(8, 5))
plt.scatter(df['MedInc'], df['MedHouseVal'], alpha=0.1, s=5, color='steelblue')
plt.axhline(5.0, color='red', linestyle='--', linewidth=1, label='Ceiling at 5.0')
plt.xlabel('Median Income ($10k)')
plt.ylabel('Median House Value ($100k)')
plt.title('Income vs. House Value')
plt.legend()
plt.tight_layout()
plt.show()

The ceiling effect is visible as a horizontal band of points at 5.0, especially at higher incomes. High-income blocks frequently hit the ceiling — the true values are hidden above it.

Geographic Analysis

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

# All blocks colored by house value
sc = axes[0].scatter(
    df['Longitude'], df['Latitude'],
    c=df['MedHouseVal'], cmap='plasma',
    alpha=0.4, s=2
)
plt.colorbar(sc, ax=axes[0], label='MedHouseVal ($100k)')
axes[0].set_title('House Value by Location')
axes[0].set_xlabel('Longitude')
axes[0].set_ylabel('Latitude')

# Same map colored by income
sc2 = axes[1].scatter(
    df['Longitude'], df['Latitude'],
    c=df['MedInc'], cmap='YlOrRd',
    alpha=0.4, s=2
)
plt.colorbar(sc2, ax=axes[1], label='MedInc ($10k)')
axes[1].set_title('Income by Location')
axes[1].set_xlabel('Longitude')
axes[1].set_ylabel('Latitude')

plt.tight_layout()
plt.show()

Info

The geographic signal is strong and non-linear. Coastal areas (Santa Barbara, Malibu, Marin County, the Peninsula south of San Francisco) cluster at high values. The Central Valley runs diagonally through the middle with low values. This pattern cannot be captured by a linear model using raw latitude and longitude — the relationship is spatial and local. Tree-based models handle this naturally.

Outlier Detection in Room Counts

# How extreme do AveRooms and AveBedrms get?
print(df[['AveRooms', 'AveBedrms', 'AveOccup']].describe())
#         AveRooms  AveBedrms   AveOccup
# count   20640.0    20640.0    20640.0
# mean        5.43       1.10       3.07
# std         2.47       0.47      10.39
# max       141.91      34.07    1243.33

# Blocks with extreme room counts
extreme_rooms = df[df['AveRooms'] > 20]
print(f"Blocks with AveRooms > 20: {len(extreme_rooms)}")
# Blocks with AveRooms > 20: 78

# What are these blocks?
print(extreme_rooms[['AveRooms', 'AveBedrms', 'AveOccup', 'Population', 'MedHouseVal']].head(5))
#      AveRooms  AveBedrms  AveOccup  Population  MedHouseVal
# ...  various extreme values

# Visualize the outliers
fig, axes = plt.subplots(1, 3, figsize=(14, 4))

for ax, col in zip(axes, ['AveRooms', 'AveBedrms', 'AveOccup']):
    ax.boxplot(df[col], vert=True, patch_artist=True,
               boxprops=dict(facecolor='steelblue', alpha=0.6))
    ax.set_title(col)
    ax.set_ylabel('Value')

plt.suptitle('Outliers in Room/Occupancy Columns', y=1.02)
plt.tight_layout()
plt.show()

Warning

A block with AveRooms = 141 is almost certainly a group-quarters situation — a dormitory, prison, or military barracks — where a small number of households have a huge shared building. These are not typical residential blocks. Leaving these extreme values in the dataset will distort linear models significantly. You will handle them in feature engineering.

Population Distribution

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

axes[0].hist(df['Population'], bins=80, color='darkorange', alpha=0.7, edgecolor='none')
axes[0].set_title('Population Distribution (Raw)')
axes[0].set_xlabel('Population')

axes[1].hist(np.log1p(df['Population']), bins=80, color='teal', alpha=0.7, edgecolor='none')
axes[1].set_title('Population Distribution (Log-transformed)')
axes[1].set_xlabel('log(1 + Population)')

plt.tight_layout()
plt.show()

print(df['Population'].skew().round(3))
# 6.175  — very right-skewed

The log transform makes Population nearly symmetric. The same applies to AveOccup.

Pairplot of Key Features

# Pairplot of the four most informative features + target
key_cols = ['MedInc', 'HouseAge', 'AveRooms', 'AveOccup', 'MedHouseVal']
sample = df[key_cols].sample(2000, random_state=42)  # sample for speed

sns.pairplot(sample, plot_kws={'alpha': 0.3, 's': 5})
plt.suptitle('Pairplot — Key Features', y=1.01)
plt.show()

EDA Summary — Key Findings

Success

What the EDA tells you before you touch a model:

1. Income dominates. MedInc alone explains ~47% of variance (r = 0.69, r² ≈ 0.47). No other single feature comes close.
2. Geography matters but is non-linear. Coastal blocks are expensive; inland blocks are cheap. Latitude and longitude carry real signal but need a non-linear model to exploit it fully.
3. The ceiling at 5.0 is a real problem. 4.7% of blocks are capped. Your model will systematically underpredict the most expensive areas.
4. Outliers exist in room counts. A small number of extreme blocks (dorms, group quarters) will pull linear models. Cap these before training.
5. Population and AveOccup are right-skewed. Log-transform them before fitting linear models.
6. No missing values. Spend your preprocessing time on feature engineering, not imputation.

---

← Back: Dataset Guide | Next: Feature Engineering →