This post is part 1 of approaching a machine learning project. Given a problem statement, we will try to see how this can be solved using machine learning and what type of algorithm is best suited to solve it. Often, we tend to jump into solving the problem without understanding the underlying complexities. This increases the chances of inaccurate predictions. This post tries to change that mindset, and helps you to approach the problem from a investigative mindset. So, let’s dive into it right away!!

This post, including the code logic has been taken from the following book, as the objective is to teach and explain whatever I learn from this book : Hands-On Machine Learning with Scikit-Learn, Keras, and TensorFlow: Concepts, Tools, and Techniques to Build Intelligent Systems 2nd Edition by  Aurélien Géron

To solve a problem, below are the important steps that one should keep in mind

1. Understand the overall context.
2. Gather relevant data.
3. Analyze and visualize the data for insights.
4. Process the data for machine learning applications.
5. Choose a model and train it.
6. Optimize your model for performance.
7. Share your solution effectively.
8.  Implement, monitor, and sustain your system.

Understand the overall context

The problem statement here is to create a model that predicts housing prices in California using census data. This dataset includes various metrics for each block group in the state, such as population, median income, and median housing price. Block groups, which typically contain 600 to 3,000 residents, represent the smallest geographical unit for which the US Census Bureau releases sample data. For simplicity, we’ll refer to these as “districts.” The model should analyze this information and accurately predict the median housing price for any district based on the other available metrics.

Framing the problem right:

• The first step is to clarify the business objective; remember, building a model isn’t necessarily the end goal. How does the company plan to utilize and benefit from this model? Understanding this is crucial, as it will influence how you define the problem, select algorithms, determine performance metrics, and allocate time for fine-tuning. You find out, that the output of your model, which predicts the median housing price for a district, will be integrated into another Machine Learning system alongside various other data points. This subsequent system will assess the investment potential of specific areas. Getting this right is vital, as it has a direct impact on revenue.
• Next, inquire about the current solution, if one exists. This can provide a benchmark for performance and valuable insights into the problem at hand. Your boss shares that district housing prices are currently estimated manually by a team of experts. They compile up-to-date information for each district, and when precise median housing prices aren’t available, they resort to complex estimation rules. This process is costly and time-consuming, and the estimates often miss the mark, sometimes by over 20%. This inefficiency is why the company believes training a model to predict median housing prices based on available data could be beneficial. The census data appears to be a valuable resource, offering median housing prices and additional metrics for thousands of districts.

With this context, you can begin designing your system. First, you need to frame the problem:

1. Is it supervised, unsupervised, or reinforcement learning?
2. Is it a classification task, a regression task, or something else?
3. Should you apply batch learning or online learning techniques?

In this case, it’s clearly a supervised learning problem since you have labeled training examples (each instance includes the expected output—the district’s median housing price). Furthermore, it qualifies as a regression task, as you are tasked with predicting a value. More specifically, this is a multiple regression problem because multiple features (such as population and median income) will inform the prediction. It is also a univariate regression problem, given that you are predicting a single value for each district; if multiple values were being predicted, it would be a multivariate regression problem. Lastly, since there isn’t a continuous stream of incoming data that requires rapid adjustments, and the dataset is manageable enough to fit in memory, a standard batch learning approach should suffice.

Selecting a Performance Measure

Root Mean Squared Error (RMSE) is a commonly used metric to evaluate the performance of regression models. It measures the average magnitude of the errors between predicted and actual values, providing insight into how well the model is performing.

The RMSE is calculated using the following steps:

• Calculate the Residuals: For each predicted value [math]\hat{Y}_i[/math] and actual value [math]Y_i[/math], compute the residual (error): [math]e_i = y_i – \hat{y}_i[/math]
• Square the Residuals: Square each of the residuals:  [math]\text{MSE} = \frac{1}{n} \sum_{i=1}^{n} e_i^2[/math], where n  is the number of observation.
• Square Root of the Mean: Finally, take the square root of the mean squared error: [math]\text{RMSE} = \sqrt{\text{MSE}} = \sqrt{\frac{1}{n} \sum_{i=1}^{n} (y_i – \hat{y}_i)^2}[/math]

Root Mean Squared Error (RMSE) gives greater weight to larger errors due to the squaring of the residuals in its calculation. Here’s a breakdown of how this works:

Squaring the Errors:
1. Residuals: For each prediction, the error is calculated as the difference between the actual value [math]Y_i[/math] and the predicted value [math]\hat{Y}_i[/math]  (i.e., [math]e_i = y_i – \hat{y}_i[/math]).
2. Squaring: Each residual is squared, which means that larger errors become disproportionately larger:
– For example, if an error is 2, squaring it gives 4, while an error of  10 results in 100. This non-linear transformation emphasizes the impact of larger errors.

Mean Calculation:
3. Mean of Squared Errors: The Mean Squared Error (MSE) is calculated by averaging these squared errors:
[math]\text{MSE} = \frac{1}{n} \sum_{i=1}^{n} e_i^2[/math]
Larger squared errors contribute significantly more to the mean compared to smaller ones.

Square Root:
4. Taking the Square Root: The final RMSE value is obtained by taking the square root of the MSE. While this step brings the RMSE back to the original units of the target variable, the influence of large errors remains due to their earlier squaring.

Summary:
Amplified Impact: Because RMSE squares the errors, larger discrepancies between predicted and actual values have a much larger impact on the final RMSE score. This means that a model with a few significant outliers can have a much higher RMSE than one with smaller, more uniformly distributed errors.
Sensitivity to Outliers: This characteristic makes RMSE particularly sensitive to outliers, making it an important consideration when evaluating model performance.

In practice, if you’re concerned about the influence of large errors or outliers, you might also consider alternative metrics like Mean Absolute Error (MAE), which treats all errors equally regardless of their size.

Mean Squared Error (MAE) :

Mean Absolute Error (MAE) is a commonly used metric for evaluating the accuracy of regression models. It measures the average magnitude of errors in a set of predictions, without considering their direction (i.e., whether they are overestimates or underestimates).

Formula:

The MAE is calculated using the following steps:

Calculate the Absolute Residuals: For each actual value [math]Y_i[/math] and the predicted value [math]\hat{Y}_i[/math], compute the absolute error:

[math]|e_i| = |y_i – \hat{y}_i|[/math]

Mean of the Absolute Errors: Take the average of these absolute errors:

[math]\text{MAE} = \frac{1}{n} \sum_{i=1}^{n} |y_i – \hat{y}_i|[/math]

where n is the number of observations.

Key Characteristics:

• Interpretability: MAE is easy to interpret since it is in the same units as the target variable. For example, if predicting prices, an MAE of $10 means the predictions are, on average, $10 off.
• Sensitivity to Outliers: Unlike Root Mean Squared Error (RMSE), MAE treats all errors equally since it does not square the residuals. This makes MAE more robust to outliers, as large errors do not disproportionately influence the overall metric.
• Linear Measurement: Because MAE uses absolute values, it provides a linear score that reflects the average error magnitude without exaggerating the influence of larger errors.

Use Cases:

MAE is particularly useful in situations where:

• You want a straightforward measure of average error.
• You want to minimize the impact of outliers.
• Interpretability is a priority, as MAE provides a clear understanding of model accuracy.

Comparison with Other Metrics:

• MAE vs. RMSE: While RMSE gives more weight to larger errors (sensitive to outliers), MAE treats all errors uniformly, which can be beneficial in specific contexts.
• MAE vs. MSE: Like RMSE, Mean Squared Error (MSE) emphasizes larger errors due to squaring, whereas MAE does not.

Summary:

Mean Absolute Error (MAE) is a valuable metric for assessing the performance of regression models, providing a clear and interpretable measure of average prediction error while being less influenced by outliers compared to metrics like RMSE.

Review the assumptions : It’s essential to identify and verify the assumptions that have been made—either by you or others—throughout the process. This step can help catch significant issues early on. For instance, your model outputs district prices that will feed into a downstream Machine Learning system. We assume these prices will be used directly, but what if the downstream system actually categorizes them into groups like “cheap,” “medium,” or “expensive”? In that case, the precise price accuracy may be irrelevant; the model only needs to assign the correct category. If that’s the situation, the problem should have been framed as a classification task rather than a regression task. You certainly don’t want to discover this after investing months in developing a regression model.

In this post, we will go over all the nitty gritty details around Multiple Linear Regression, but first let’s start with understanding what is Multiple Linear Regression.

Multiple linear regression is a statistical technique used to model the relationship between multiple independent variables (also called predictors or features) and a dependent variable (also called the outcome or response variable). It assumes a linear relationship between the predictors and the outcome variable.

Model Representation

The multiple linear regression model can be represented as:

[math]Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \ldots + \beta_n X_n + \varepsilon[/math]

Where:

• [math]Y[/math] is the dependent variable (outcome).
• [math]X_1, X_2, \ldots, X_n[/math]are the independent variables (predictors).
• [math]\beta_0[/math] is the intercept (constant term).
• [math]\beta_1, \beta_2, \ldots, \beta_n[/math] are the coefficients (parameters) representing the change in [math]Y[/math] for a unit change in each corresponding [math]X[/math].
• [math]\varepsilon[/math] is the error term, representing the difference between the observed and predicted values.

Mathematical Derivation

The objective of Multiple Linear Regression is to estimate the coefficients [math]\beta_1, \beta_2, \ldots, \beta_n[/math] that minimize the sum of squared differences between the observed and predicted values. This is typically done using the method of ordinary least squares (OLS).

The OLS estimates for the coefficients can be obtained by minimizing the sum of squared residuals (errors):

[math]\sum_{i=1}^{n} (Y_i – (\beta_0 + \beta_1 X_{i1} + \beta_2 X_{i2} + \ldots + \beta_n X_{in}))^2[/math]

To minimize this expression, the partial derivatives of the sum of squared residuals with respect to each coefficient are equated to zero:

[math]\frac{\partial}{\partial \beta_0} \sum_{i=1}^{n} (Y_i – (\beta_0 + \beta_1 X_{i1} + \beta_2 X_{i2} + \ldots + \beta_n X_{in}))^2 = 0[/math]

[math]\frac{\partial}{\partial \beta_1} \sum_{i=1}^{n} (Y_i – (\beta_0 + \beta_1 X_{i1} + \beta_2 X_{i2} + \ldots + \beta_n X_{in}))^2 = 0[/math]

[math]\vdots[/math]

[math]\frac{\partial}{\partial \beta_n} \sum_{i=1}^{n} (Y_i – (\beta_0 + \beta_1 X_{i1} + \beta_2 X_{i2} + \ldots + \beta_n X_{in}))^2 = 0[/math]

Before going any further, we should understand why did we do partial derivative on the [math]\sum_{i=1}^{n} (Y_i – (\beta_0 + \beta_1 X_{i1} + \beta_2 X_{i2} + \ldots + \beta_n X_{in}))^2[/math] expression and equate it to zero

As we already saw in the above section [math]\sum_{i=1}^{n} (Y_i – (\beta_0 + \beta_1 X_{i1} + \beta_2 X_{i2} + \ldots + \beta_n X_{in}))^2[/math] expression represents the residual sum of squares (RSS) in the context of linear regression, which quantifies the discrepancy between the observed values [math]Y_i[/math] and the predicted values [math]\hat{Y}_i[/math] generated by the regression model

RSS = [math]\sum_{i=1}^{n} (Y_i – (\beta_0 + \beta_1 X_{i1} + \beta_2 X_{i2} + \ldots + \beta_n X_{in}))^2[/math]

Or it can also be represented as

RSS = [math]\sum_{i=1}^{n} (Y_i – X_i \beta)^2[/math]

The goal in linear regression is to find the values of the coefficients [math]\beta_1, \beta_2, \ldots, \beta_n[/math]  that minimize the RSS. This is typically done by finding the values of the coefficients that make the partial derivatives of the RSS with respect to each coefficient equal to zero as show in the equations above. These values are then the critical points where the RSS might be at a minimum.

Now, let’s understand how the partial derivative works:

• Partial Derivatives: Taking the partial derivatives of the RSS with respect to each coefficient [math]\beta_j \quad \text{where} \quad j = 0, 1, \ldots, n[/math] provides information about how the RSS changes with respect to changes in each coefficient individually. Setting these partial derivatives to zero helps identify candidate solutions where the RSS might be minimized. That’s the reason why we are taking the partial derivative of RSS with respect to each coefficient and equating it to 0

[math]\frac{\partial}{\partial \beta_0} \sum_{i=1}^{n} (Y_i – (\beta_0 + \beta_1 X_{i1} + \beta_2 X_{i2} + \ldots + \beta_n X_{in}))^2 = 0[/math]

[math]\frac{\partial}{\partial \beta_1} \sum_{i=1}^{n} (Y_i – (\beta_0 + \beta_1 X_{i1} + \beta_2 X_{i2} + \ldots + \beta_n X_{in}))^2 = 0[/math]

[math]\vdots[/math]

[math]\frac{\partial}{\partial \beta_n} \sum_{i=1}^{n} (Y_i – (\beta_0 + \beta_1 X_{i1} + \beta_2 X_{i2} + \ldots + \beta_n X_{in}))^2 = 0[/math]

Vector Representation

Now that we have understood the mathematical model representation, lets try to interpret it in the form of vector representation, so that it will be easier for us to translate it into code while modeling the data

Consider the Multiple regression model : [math]Y = X\beta + \epsilon[/math]

RSS = [math]\sum_{i=1}^{n} (Y_i – X_i \beta)^2[/math]

To represent RSS in matrix form, let’s denote the following:

• [math]Y[/math] as the vector of observed values of the dependent variable (of length n).
• [math]X[/math] as the matrix of predictor variables (of size n \times (p+1), where p is the number of predictors, and the 1 represents the first column of X whichis a vector of ones for the intercept term).
• [math]\hat{Y}[/math] as the vector of predicted values of the dependent variable.

Then, the predicted value [math]\hat{Y}[/math] can be expressed as [math]\hat{Y}[/math]  = [math]X \beta[/math], where [math]X \beta[/math] represents the matrix multiplication of [math]X[/math] and [math]\beta[/math]

Now we could define a residual vector [math]\varepsilon[/math] which denotes the vector difference between observed values of [math]Y[/math] and the predicted values of [math]\hat{Y}[/math], i.e., [math]\varepsilon[/math] = [math]Y[/math]  – [math]\hat{Y}[/math]

Using the above equation we can rewrite the equation of RSS as the sum of square of the residual vector [math]\varepsilon[/math]

[math]RSS = \sum_{i=1}^{n} \epsilon_i^2[/math]

In matrix form, this can be written as the dot product of the transpose math]\varepsilon[/math] and [math]\varepsilon[/math]:

[math]RSS = \epsilon^T \epsilon[/math]

Expanding this expression using the definitions of [math]\varepsilon[/math] and [math]\hat{Y}[/math], we get:

[math]RSS = (Y – X\beta)^T (Y – X\beta)[/math]

Expanding and simplifying the above expression:

[math]RSS = Y^T Y – 2\beta^T X^T Y + \beta^T X^T X \beta[/math]

To find the coefficients [math]\beta[/math] that minimize the RSS, we take the derivative of RSS with respect to [math]/beta[/math],set it equal to zero, and solve for [math]/beta[/math]. This is the same method we used above in our partial derivative section. The only difference was we are solving this in a vectorized format and the earlier version was in scaler format.

[math]\frac{\partial RSS}{\partial \beta} = -2X^T Y + 2X^T X\beta = 0[/math]

Now, solving the above equation for [math]\beta[/math] gives us:

[math](X^T X) \beta = X^T Y[/math]
[math]\hat{\beta} = (X^T X)^{-1} X^T y[/math]

where [math]\hat{\beta}[/math] is the vector of estimated coefficients.

Now that we have understood the Model representation and the mathematical derivation of Multiple Linear Regression, lets implement a model in code. I will first write each block of code and provide explaination and later consolidate the whole code

import numpy as np
class MultipleLinearRegression:
    def __init__(self):
        self.intercept_ = None
        self.coef_ = None

This above block of code defines a Python class named MultipleLinearRegression. It initializes two instance variables intercept_ and coef_ to None.

def fit(self, X, y):
        ones = np.ones((X.shape[0], 1))
# Adding a column of ones for the intercept term
        X = np.concatenate((ones, X), axis=1)  
# Compute coefficients using closed-form solution (OLS)
        self.coefficients = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(y)
        self.intercept_ = self.coefficients[0]
        self.coef_ = self.coefficients[1:]

The above method fit is used to train the Multiple Linear Regression model. It takes two arguments X (the input features) and y (the target variable). It first creates a column vector of ones using np.ones((X.shape[0], 1)). This vector is added as the first column to the input features X using np.concatenate. This is done to account for the intercept term in the linear regression model.

Next, it computes the coefficients of the linear regression model using the closed-form solution, also known as Ordinary Least Squares (OLS). It calculates the coefficients by first multiplying the transpose of X with X representing (X.T.dot(X)), then taking the inverse of the resulting matrix (np.linalg.inv) , and finally multiplying the result with the transpose of X and y to get the coefficients. This is something which we already did during the vector representation section.
The intercept term is stored in self.intercept_, and the remaining coefficients are stored in self.coef_.

    def predict(self, X):
        ones = np.ones((X.shape[0], 1))
        X = np.concatenate((ones, X), axis=1)  # Adding a column of ones for the intercept term
        return X.dot(self.coefficients)

This method predict is used to make predictions using the trained linear regression model.
Similar to the fit method, it adds a column of ones to the input features to account for the intercept term.
It then calculates the predicted values (y_pred) using the dot product of the input features X and the coefficients self.coefficients

# Example usage:
# Generate some sample data
np.random.seed(0)
X = np.random.rand(100, 2)  # 100 samples, 2 features
# True relationship: y = 2*x1 + 3*x2 + noise
y = 2*X[:,0] + 3*X[:,1] + np.random.randn(100)  
# Split the data into train and test sets
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)
# Fit the multiple linear regression model
model = MultipleLinearRegression()
model.fit(X_train, y_train)
# Make predictions on the test set
y_pred = model.predict(X_test)
# Evaluate the model
from sklearn.metrics import r2_score
accuracy = r2_score(y_test, y_pred)
print('Accuracy (R^2 score):', accuracy)
print('Intercept:', model.intercept_)
print('Coefficients:', model.coef_)

This block of code demonstrates the usage of the MultipleLinearRegression class.It generates some sample data X and y where the relationship between features and the target variable y is known.
Then, it splits the data into training and testing sets using train_test_split from scikit-learn.It fits the MultipleLinearRegression model to the training data using the fit method.After that, it makes predictions on the test set using the predict method.Finally, it evaluates the model’s performance using the R^2 score and prints the result.

Finally running the whole code

import numpy as np
class MultipleLinearRegression:
    def __init__(self):
        self.intercept_ = None
        self.coef_ = None
    def fit(self, X, y):
        ones = np.ones((X.shape[0], 1))
        X = np.concatenate((ones, X), axis=1)  # Adding a column of ones for the intercept term
        # Compute coefficients using closed-form solution (OLS)
        self.coefficients = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(y)
        self.intercept_ = self.coefficients[0]
        self.coef_ = self.coefficients[1:]
    def predict(self, X):
        ones = np.ones((X.shape[0], 1))
        X = np.concatenate((ones, X), axis=1)  # Adding a column of ones for the intercept term
        return X.dot(self.coefficients)
# Example usage:
# Generate some sample data
np.random.seed(0)
X = np.random.rand(100, 2)  # 100 samples, 2 features
y = 2*X[:,0] + 3*X[:,1] + np.random.randn(100)  # True relationship: y = 2*x1 + 3*x2 + noise
# Split the data into train and test sets
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)
# Fit the multiple linear regression model
model = MultipleLinearRegression()
model.fit(X_train, y_train)
# Make predictions on the test set
y_pred = model.predict(X_test)
# Evaluate the model
from sklearn.metrics import r2_score
accuracy = r2_score(y_test, y_pred)
print('Accuracy (R^2 score):', accuracy)
print('Intercept:', model.intercept_)
print('Coefficients:', model.coef_)
# Plotting the model's predictions vs. actual values
plt.scatter(y_test, y_pred, color='blue')
plt.plot(y_test, y_test, color='red', linestyle='--')  # Plotting the y = x line for reference
plt.title('Model Predictions vs. Actual Values')
plt.xlabel('Actual Values')
plt.ylabel('Predicted Values')
plt.show()
Accuracy (R^2 score): 0.43228286040209773
Intercept: 0.04493496059558341
Coefficients: [1.56277657 3.09631373]
​

References:

• Regression Analysis
• Regression Equation