Linear Regression is one of the simplest and most important machine learning algorithms for predicting continuous numerical values.

In this project, I used the California Housing dataset available in Scikit-Learn to build a Linear Regression model that predicts the median house value of a district.

The goal was not only to train the model, but also to understand what happens at each stage:

1. Load the dataset
2. Explore the data
3. Check correlations
4. Split the data into training and testing sets
5. Train Linear Regression
6. Understand the learned coefficients
7. Make predictions
8. Evaluate the model
9. Analyze residuals
10. Investigate how errors change across prediction ranges

1. What is Linear Regression?

Linear Regression is a supervised machine learning algorithm used to predict a continuous numerical value.

For example:

• Predict house price
• Predict salary
• Predict temperature
• Predict sales
• Predict electricity consumption

The basic idea is to find a mathematical relationship between input features and the target.

For one feature:

$\hat{y} = \beta_0 + \beta_1x$

Where:

• $\hat{y}$ = predicted value
• $x$ = input feature
• $\beta_0$ = intercept
• $\beta_1$ = coefficient

With multiple features:

$\beta_0+\beta_1x_1+\beta_2x_2+\cdots+\beta_nx_n$

This is called Multiple Linear Regression.

2. Project Goal

For this project, the goal is:

Predict the median house value using demographic and geographic information about California districts.

We will use the California Housing dataset provided by Scikit-Learn.

3. Loading the Dataset

First, import the required libraries.

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

from sklearn.datasets import fetch_california_housing

Load the dataset:

housing = fetch_california_housing()

The dataset contains:

• housing.data → input features
• housing.feature_names → feature names
• housing.target → target values

Convert the features into a Pandas DataFrame:

X = pd.DataFrame(
    housing.data,
    columns=housing.feature_names
)

Convert the target into a Pandas Series:

y = pd.Series(
    housing.target,
    name="Price"
)

Now:

X.shape

gives:

(20640, 8)

And:

y.shape

gives:

(20640,)

So we have:

• 20,640 observations
• 8 input features
• 1 target variable

4. Understanding the Features

The dataset contains the following features:

The target is:

Price

The target represents the median house value in units of $100,000.

Therefore:

Price = 2.0

means approximately:

$200,000

Similarly:

Price = 5.0

means approximately:

$500,000

5. Exploratory Data Analysis

Before training a machine learning model, we should understand the data.

A useful first step is:

X.describe()

This gives statistics such as:

• count
• mean
• standard deviation
• minimum
• maximum
• quartiles

For example, MedInc has approximately:

Mean = 3.87
Minimum = 0.50
Maximum = 15.00

Some variables also contain extreme values.

For example:

AveRooms maximum ≈ 141.91
AveOccup maximum ≈ 1243.33

These are unusually large compared with their average values and are worth investigating because extreme observations can influence Linear Regression.

6. Checking Missing Values

Before training the model, we should check whether the features contain missing values.

X.isnull().sum()

The California Housing dataset contains no missing values in these features.

This means we do not need to perform missing-value imputation for this baseline model.

7. Correlation Analysis

Correlation helps us understand the linear relationship between two variables.

We can calculate the correlation between every feature and the target:

correlation = X.copy()

correlation["Price"] = y

correlation.corr()["Price"].sort_values(
    ascending=False
)

The results were:

8. How to Interpret Correlation

The correlation coefficient generally ranges from:

$-1 \leq r \leq 1$

Positive correlation

If:

$r > 0$

then the two variables tend to increase together.

For example:

MedInc → Price
r = 0.688

This indicates a relatively strong positive linear relationship.

In simple terms:

Districts with higher median income tend to have higher median house values.

Negative correlation

If:

$r < 0$

then one variable tends to increase as the other decreases.

For example:

Latitude → Price
r = -0.144

This is a weak negative linear relationship.

Correlation close to zero

If:

$r \approx 0$

there is little linear relationship between the variables individually.

For example:

Population → Price
r ≈ -0.025

This does not necessarily mean Population is useless.

A feature may have a weak individual correlation but still provide useful information when combined with other features.

Also:

Correlation does not imply causation.

9. Splitting the Dataset

We should not train and evaluate the model using exactly the same data.

Why?

Because we want to know whether the model can make predictions on unseen data.

We split the dataset into:

• Training data
• Testing data

Using:

from sklearn.model_selection import train_test_split

X_train, X_test, y_train, y_test = train_test_split(
    X,
    y,
    test_size=0.2,
    random_state=42
)

What does test_size=0.2 mean?

It means:

80% → Training
20% → Testing

Our dataset contains 20,640 observations.

Therefore:

Training = 16,512
Testing  = 4,128

The resulting shapes were:

Original dataset:   (20640, 8)
Training features:  (16512, 8)
Testing features:   (4128, 8)

Training target:    (16512,)
Testing target:     (4128,)

10. Why Do We Use random_state=42?

We used:

random_state=42

This controls the random splitting process.

Without it, every execution could produce a different train/test split.

With:

random_state=42

we get the same split every time.

The number 42 itself has no special mathematical meaning here. It is simply a commonly used fixed seed.

11. Training the Linear Regression Model

Now we can create our Linear Regression model.

from sklearn.linear_model import LinearRegression

model = LinearRegression()

Train it:

model.fit(X_train, y_train)

The .fit() method learns the coefficients from the training data.

Conceptually, the model is trying to find:

$\beta_0+\beta_1x_1+\beta_2x_2+\cdots+\beta_8x_8$

so that its predictions are as close as possible to the actual target values.

12. Understanding the Intercept

We obtained:

Intercept = -37.0233

The intercept is:

$\beta_0 = -37.0233$

It is the predicted target when all input features are zero.

Our model therefore starts with:

$\hat{y} = -37.0233 + …$

In real-world interpretation, the intercept does not necessarily have a meaningful physical interpretation because a situation where all eight features are zero may not represent a realistic California district.

13. Understanding the Coefficients

The learned coefficients were:

The learned equation is approximately:

$\hat y = -37.0233 +0.4487(MedInc) +0.0097(HouseAge) -0.1233(AveRooms) +0.7831(AveBedrms) -0.00000203(Population) -0.0035(AveOccup) -0.4198(Latitude) -0.4337(Longitude)$

14. How to Interpret a Coefficient

Consider:

MedInc coefficient = 0.4487

This means:

Holding the other features constant, a one-unit increase in MedInc is associated with an increase of approximately 0.4487 in the predicted target.

Because the target is measured in $100,000:

0.4487 \times 100{,}000 \approx \44{,}867$

So, approximately:

A one-unit increase in MedInc, while keeping the other model features fixed, corresponds to about a $44,867 increase in predicted house value.

The phrase "holding other features constant" is important because this is a multiple linear regression model.

15. Why Can a Coefficient Be Negative Even When Correlation Is Positive?

This is an important concept.

We saw:

AveRooms correlation with Price = +0.152

but the Linear Regression coefficient was:

AveRooms coefficient = -0.123

At first, this may seem contradictory.

But correlation looks at the relationship between:

AveRooms ↔ Price

individually.

The regression coefficient looks at:

AveRooms ↔ Price

while accounting for the other features in the model.

The input features themselves can be correlated with each other.

For example:

MedInc
   ↕
AveRooms
   ↕
AveBedrms

Because several predictors contain overlapping information, the regression coefficients can change direction after the other variables are included.

This is one reason why we later want to investigate multicollinearity.

16. Making Predictions

After training the model, we use the test data to make predictions.

y_pred = model.predict(X_test)

Here:

X_test

contains houses the model has never seen during training.

The model produces:

y_pred

which contains the predicted prices.

17. Comparing Actual and Predicted Values

We can create a DataFrame:

comparison = pd.DataFrame({
    "Actual Price": y_test.values,
    "Predicted Price": y_pred
})

The first few predictions looked like:

Because the target is in units of $100,000, we can convert the values to dollars:

comparison["Actual Price ($)"] = (
    comparison["Actual Price"] * 100000
)

comparison["Predicted Price ($)"] = (
    comparison["Predicted Price"] * 100000
)

For example:

Actual = 5.0

means:

$500,000

while:

Predicted = 2.7097

means approximately:

$270,966

So the model significantly underpredicted that particular house value.

18. What is a Residual?

A residual is the difference between the actual value and the predicted value.

$Residual = Actual - Predicted$

or:

$e_i = y_i - \hat{y}_i$

Example

Suppose:

Actual = 3.0
Predicted = 2.5

Then:

$Residual = 3.0 - 2.5 = 0.5$

The model underpredicted.

If:

Actual = 2.0
Predicted = 2.5

then:

$Residual = 2.0 - 2.5 = -0.5$

The model overpredicted.

19. Evaluating the Model

We used four common regression metrics:

from sklearn.metrics import (
    mean_absolute_error,
    mean_squared_error,
    r2_score
)

import numpy as np

y_pred = model.predict(X_test)

mae = mean_absolute_error(y_test, y_pred)
mse = mean_squared_error(y_test, y_pred)
rmse = np.sqrt(mse)
r2 = r2_score(y_test, y_pred)

print("Model Evaluation")
print("----------------")
print("MAE :", mae)
print("MSE :", mse)
print("RMSE:", rmse)
print("R²  :", r2)

print("\nIn dollars:")
print("MAE :", mae * 100000)
print("RMSE:", rmse * 100000)

Our results:

MAE  = 0.5332
MSE  = 0.5559
RMSE = 0.7456
R²   = 0.5758

In dollars:

MAE  ≈ $53,320
RMSE ≈ $74,558

20. Mean Absolute Error — MAE

MAE measures the average absolute prediction error.

$MAE = \frac{1}{n} \sum_{i=1}^{n} |y_i-\hat{y}_i|$

Our result:

$MAE = 0.5332$

Since the target is measured in $100,000:

0.5332 \times 100{,}000 \approx \53{,}320$

So:

On average, the model's predictions differ from the actual house values by about $53,320 in absolute terms.

MAE is easy to interpret because it is in the same unit as the target.

21. Mean Squared Error — MSE

MSE calculates the average squared error.

$MSE = \frac{1}{n} \sum_{i=1}^{n} (y_i-\hat{y}_i)^2$

Our result:

$MSE = 0.5559$

The important point is that errors are squared.

Therefore, large errors receive much more weight.

For example:

$1^2=1$

but:

$5^2=25$

So a large prediction error can have a much bigger effect on MSE.

22. Root Mean Squared Error — RMSE

RMSE is simply the square root of MSE.

$RMSE = \sqrt{MSE}$

Therefore:

$RMSE = \sqrt{0.5559} \approx 0.7456$

In dollars:

0.7456 \times 100{,}000 \approx \74{,}558$

So our RMSE is approximately:

$74,558

RMSE is particularly useful when we want large errors to receive more attention.

23. Why Is RMSE Larger Than MAE?

Our metrics were:

MAE  ≈ 0.533
RMSE ≈ 0.746

RMSE is larger because it squares the errors before averaging them.

Therefore, unusually large errors affect RMSE strongly.

We later found one residual as large as:

−9.875

which corresponds to an error of approximately:

9.875 \times 100{,}000 \approx \987{,}533$

Such a large error can significantly increase RMSE.

24. R² Score

The R² score tells us how much of the variation in the target is explained by the model.

Our result:

$R^2 = 0.5758$

or approximately:

57.6%

So we can say:

The Linear Regression model explains approximately 57.6% of the variation in median house values on this test set.

Important:

It does not mean:

"The model predicts 57.6% of houses correctly."

This is a regression problem, not a classification problem.

25. Summary of Model Performance

This gives us a reasonable baseline model.

The important thing is not simply whether the number is "good" or "bad."

The baseline gives us something that we can later compare against more advanced models.

26. Residual Analysis

After calculating metrics, we should investigate how the model is making its mistakes.

First calculate the residuals:

residuals = y_test - y_pred

Then plot residuals against predicted values:

plt.figure(figsize=(8, 5))

plt.scatter(
    y_pred,
    residuals,
    alpha=0.4
)

plt.axhline(
    0,
    linestyle="--"
)

plt.xlabel("Predicted Price")
plt.ylabel("Residual")
plt.title("Residuals vs Predicted Price")

plt.show()

27. What Should a Good Residual Plot Look Like?

Ideally, residuals should be:

• centered around zero
• randomly scattered
• without obvious curves or patterns
• without the spread changing dramatically across predictions

Conceptually:

Residual
   |
 + |   •    •      •
   |      •    •
 0 |-------------------------
   |  •     •    •     •
 - |      •      •
   |
   +------------------------> Predicted Price

A clear pattern can indicate that Linear Regression is missing some structure in the data.

28. Mean Residual

We calculated:

residuals = y_test - y_pred

print("Mean residual:", residuals.mean())

Our result was:

Mean residual: 0.003479

This is extremely close to zero.

That is encouraging because positive and negative residuals roughly balance each other.

In other words:

The model does not show a strong overall tendency to systematically overpredict or underpredict.

However, a mean residual close to zero does not mean that individual predictions are accurate.

Large positive and negative errors can cancel each other.

29. Minimum and Maximum Residual

We calculated:

print("Minimum residual:", residuals.min())
print("Maximum residual:", residuals.max())

Results:

Minimum residual = -9.8753
Maximum residual =  4.1484

Because:

$Residual=Actual−Predicted$

a negative residual means the model overpredicted.

Therefore:

-9.8753

represents a very large overprediction.

In dollars:

9.8753 \times 100{,}000 \approx \987{,}530$

The maximum residual:

4.1484

represents an underprediction of approximately:

4.1484 \times 100{,}000 \approx \414{,}840$

So the model has some very large individual prediction errors.

30. Largest Absolute Residuals

We also calculated the observations with the largest absolute errors:

print(
    residuals.abs()
    .sort_values(ascending=False)
    .head(10)
)

The result was:

1979     9.875331
6688     4.148388
10574    3.884783
12389    3.680056
19542    3.632784
4548     3.561215
459      3.452480
15652    3.357140
12069    3.331976
4630     3.288854

The largest absolute residual was:

9.8753

or approximately:

$987,533

This is an extremely large individual prediction error.

This also helps explain why our RMSE is substantially larger than our MAE.

31. Correlation Between Predictions and Residuals

We calculated:

print(
    "Correlation between predicted values and residuals:",
    np.corrcoef(y_pred, residuals)[0, 1]
)

Result:

-0.06205

This is very close to zero.

Therefore:

There is little evidence of a strong linear relationship between predicted values and residuals.

This is generally encouraging.

However, correlation alone cannot detect every possible pattern.

For example, residuals could follow a curved relationship while still having a correlation close to zero.

That's why the residual plot remains important.

32. Does the Model Make Larger Errors for Expensive Houses?

We wanted to investigate another important question:

Does prediction error increase as the predicted house value increases?

We created a residual analysis DataFrame:

residual_analysis = pd.DataFrame({
    "Actual": y_test.values,
    "Predicted": y_pred,
    "Residual": residuals
})

residual_analysis["Absolute_Error"] = (
    residual_analysis["Residual"].abs()
)

Then divided predictions into four groups:

residual_analysis["Prediction_Group"] = pd.qcut(
    residual_analysis["Predicted"],
    q=4
)

Finally:

residual_analysis.groupby(
    "Prediction_Group",
    observed=True
)["Absolute_Error"].agg(
    ["mean", "median", "max"]
)

Our results:

33. Interpreting Error by Prediction Range

Remember that the target is measured in $100,000.

Therefore, the mean absolute errors are approximately:

There is a clear pattern:

$Error \uparrow \quad\text{as}\quad Predicted\ Price \uparrow$

In simple words:

The model tends to make larger absolute errors when predicting higher house values.

34. What Is Heteroscedasticity?

When the variability of residuals changes across the range of predictions, we may have heteroscedasticity.

For example:

Homoscedasticity

Low price              High price

small errors            small errors
   • • •                   • • •
   • • •                   • • •

The error spread stays roughly constant.

Heteroscedasticity

Low price              High price

small errors            large errors
   • •                     •
   • •                  •     •
   • •               •     •       •

The error spread becomes larger.

Our group analysis suggests that the model may exhibit this behavior because:

Mean absolute error:

0.385
  ↓
0.466
  ↓
0.593
  ↓
0.690

The median error also increases:

0.308
  ↓
0.378
  ↓
0.479
  ↓
0.552

This is important because the increasing median error suggests that the pattern is not caused only by one extreme outlier.

35. What Have We Learned So Far?

Our Linear Regression baseline gives us the following picture.

What is working?

1. The model explains a substantial amount of variation.

$R^2 \approx 0.576$

2. The average residual is close to zero.

$Mean\ Residual \approx 0$

3. Prediction/residual correlation is close to zero.

$r \approx -0.062$

What needs investigation?

1. Some predictions have extremely large errors.

Largest absolute residual:

9.875 \approx \987{,}533$

2. Error increases for higher predicted prices.

Mean absolute error increased from:

$0.385 \rightarrow 0.690$

3. Some regression coefficients are surprising.

For example:

AveRooms correlation with Price = +0.152
AveRooms regression coefficient = -0.123

This suggests that relationships among the input features may be influencing the learned coefficients.

36. Current Baseline

At this point, our baseline Linear Regression model can be summarized as:

Dataset
   ↓
California Housing
   ↓
20,640 observations
   ↓
8 features
   ↓
80/20 train-test split
   ↓
Linear Regression
   ↓
Predictions
   ↓
Evaluation
   ↓
Residual Analysis

Baseline performance

MAE  ≈ $53,320
RMSE ≈ $74,558
R²   ≈ 0.576

The model is useful as a baseline because we now have a reference point.

Any improved model should be evaluated using the same test data and comparable metrics.

37. What Comes Next?

The next logical step is to investigate multicollinearity.

Multicollinearity occurs when input features contain highly overlapping information.

This can make Linear Regression coefficients unstable or difficult to interpret.

For example, we saw that:

AveRooms

had a positive raw correlation with house price but a negative regression coefficient.

One possible reason is that AveRooms is related to other predictors.

We can investigate this using Variance Inflation Factor (VIF).

The next step is therefore:

Calculate VIF for all eight features and determine whether multicollinearity is affecting our Linear Regression model.

After that, we can investigate whether feature transformations or regularized models such as Ridge and Lasso Regression improve the baseline.