import pandas as pd
# Importing eshop_revenues
eshop_revenues = pd.read_csv("https://raw.githubusercontent.com/Datakortex/Datasets/refs/heads/main/eshop_revenues.csv")Machine Learning with Linear Regression
Introduction
In Chapters Correlation and Simple Linear Regression and Multiple Linear Regression and Statistical Inference, we discussed how linear regression can be used for statistical inference. Specifically, we showed that, under the Gauss–Markov assumptions, the estimators are BLUE (Best Linear Unbiased Estimators), and we can use them to assess whether the association between an independent variable and the dependent variable is statistically significant.
However, linear regression is not limited to statistical inference; it can also be applied as a machine learning tool. In this chapter, we focus on using linear regression from a predictive perspective. We demonstrate how to train a model on one dataset and evaluate its performance on a new, unseen dataset. Along the way, we illustrate that the model optimized for prediction may differ from the one used for statistical inference, depending on the goal, and how choices such as adding predictors or interactions can impact predictive performance.
This chapter also introduces scikit-learn, the most widely used machine learning library in Python. While the underlying linear regression model remains the same, scikit-learn provides a different workflow that emphasizes model training, prediction, and performance evaluation. This workflow forms the foundation for many of the machine learning methods presented in subsequent chapters.
Simple Linear Regression in Machine Learning
Let’s start by loading Pandas and importing the eshop_revenues dataset:
Technically, we should first proceed with exploratory data analysis, data pre-processing and feature engineering. We are already familiar with this dataset though from previous chapters and, to keep things simple, we focus immediately on the training-test phase. Since our task is to make good predictions, we need to make sure that these predictions will occur on a dataset that our model has not seen in the estimation phase (training).
Before we use the train_test_split() function though, we separate the predictors from the target variable by creating two different objects: X for the predictors and y for the target variable. We can do this by creating Pandas Series/DataFrames from the relevant columns of our dataset. Then, we use the train_test_split() function, which takes the predictors (X) and the target variable (y) as inputs and randomly partitions them into training and test sets. In the code below, the argument test_size = 0.20 specifies that 20% of the observations should be allocated to the test set, while the remaining 80% are used for training the model. Also, the argument random_state = 42 sets the seed of the random number generator, ensuring that the same split is obtained each time the code is executed.
from sklearn.model_selection import train_test_split
# Separating the predictors from the target variable
X = eshop_revenues[["Ad_Spend"]]
y = eshop_revenues["Revenue"]
# Splitting data into training and test sets
X_train, X_test, y_train, y_test = train_test_split(
X, y, test_size = 0.20, random_state = 42
)The function train_test_split() returns four objects:
X_trainandy_train, which contain the predictor and target values used to train the model.X_testandy_test, which contain the corresponding data reserved for model evaluation.
By keeping the test observations separate during model training, we can obtain a more realistic assessment of how well the model is likely to perform on new, unseen data.
Note that because scikit-learn expects the predictors to be provided as a two-dimensional structure—where rows represent observations and columns represent predictor variables—we use double brackets when selecting Ad_Spend. This keeps X as a DataFrame with one column instead of converting it into a Pandas Series. The target variable y, on the other hand, is stored as a Pandas Series because it represents a single outcome that we want the model to predict.
We can inspect these four different datasets using the head() method:
# Printing the first few rows of the training predictors
X_train.head()| Ad_Spend | |
|---|---|
| 70 | 77.22 |
| 78 | 222.29 |
| 47 | 19.87 |
| 0 | 114.19 |
| 12 | 68.75 |
Notice how the index values are randomly split between the training and test sets. The index values are not rearranged or reset; they continue to identify the original observations from the dataset. Also, notice how X_train contains only the values of the predictor variable Ad_Spend. The corresponding target values are stored separately in y_train. The index acts as a link between the two objects: a row with the same index in X_train and y_train refers to the same original observation.
# Printing the first few rows of the training target variable
y_train.head()70 663.75
78 14802.05
47 839.70
0 29089.69
12 15667.10
Name: Revenue, dtype: float64Now, we want to fit a linear regression model on the training set, using X_train and y_train as inputs. Before doing this, we need to import the LinearRegression() class from sklearn.linear_model module. We first create an object representing the model and then use the fit() method to estimate the model parameters from the training data.
from sklearn.linear_model import LinearRegression
# Creating a linear regression model
lm_model_simple = LinearRegression()
# Estimating the model parameters using the training data
lm_model_simple.fit(X_train, y_train)The slope and intercept of the estimated linear regression model can be accessed using the coef_ and intercept_ attributes:
# Printing the slope coefficient
lm_model_simple.coef_array([50.32577847])# Printing the intercept
lm_model_simple.intercept_np.float64(6323.339015187221)These estimates differ from those obtained using the full dataset. Scikit-learn’s LinearRegression() does not calculate or store p-values for coefficients. It gives us the estimated parameters (slope and intercept), but not the statistical inference output. Had we used the ols() function from the statsmodels package though, we would find that both the intercept and the slope are statistically significantly different from 0 at the 1% significance level, suggesting evidence of a linear association between the predictor and the response variable. However, our focus here is not hypothesis testing but prediction accuracy. Mathematically, we estimated the following:
\[\widehat{\text{Revenue}} = 6323.34 + 50.33 \times \text{AdSpend}\]
To test the model’s accuracy, we use it to make predictions on the test set. This is possible because the test set contains all the variables (features) from the original dataset. For example, if a test set observation has Ad_Spend = 100, the predicted revenue will be:
\[6323.34 + 50.33 \times 100 = 11356.34\]
As with every regression task, we will never predict the actual revenues exactly; the goal is to produce predictions as close as possible to the true values.
We can use the lm_model_simple object to make predictions on the test set with the predict() method. This method takes the predictor values as input and returns the estimated values of the target variable. In our case, we provide X_test, which contains the Ad_Spend values that were not used during model training. The model then uses the estimated slope and intercept to calculate the predicted revenue values.
# Generating predictions for the test set
y_pred = lm_model_simple.predict(X_test)
# Printing the first few values
y_pred[:5]array([16446.87261144, 18814.70048828, 21741.64776387,
13050.88908055, 9798.83727606])The resulting values in y_pred represent the model’s predicted revenues for the observations in X_test. Because X_test and y_test have the same index values, we can compare these predictions with the actual observed revenues stored in y_test to evaluate the performance of the model.
At this point, the key question is how we can make such an assessment. We discuss more advanced evaluation metrics in later chapters, but one simple and intuitive metric is the mean absolute error (MAE). This is the average absolute difference between the predicted revenues and the actual revenues:
\[\text{MAE} = \frac{1}{N} \sum^N_{i = 1}|y_i - \hat{y}_i|\]
This metric is expressed in the same units as the target variable (revenue) and measures, on average, how far predictions are from the true values. Therefore, the lower the value of MAE, the better our predictions are on average. Because it uses absolute values, over-predictions (predicted revenues higher than actual) and under-predictions (predicted revenues lower than actual) do not cancel each other out, which makes the metric easy to interpret.
To evaluate the predictive performance of our model, we can calculate the MAE on the test set using scikit-learn. This is done by comparing the true values from y_test with the predicted values generated by our model for the test predictors X_test. In other words, we first generate predictions using the fitted model and then compare those predictions to the actual observed outcomes.
Before doing so, we import the mean_absolute_error() function from sklearn.metrics. This function implements the MAE formula and allows us to compute it directly from the true and predicted values:
from sklearn.metrics import mean_absolute_error
# Calculating the mean absolute error
mae = mean_absolute_error(y_test, y_pred)
# Printing the result
mae9699.194593476483The result shows that, on average, our predictions differ from the actual revenues by approximately €9,700. Whether this is a “good” result depends on the context: the scale of revenues, the business application, and the cost of prediction errors. In machine learning, model performance is always relative to the problem being solved.
As a baseline comparison, we can consider a very simple model: using only the average revenue of the training set as the prediction for every observation in the test set. In this case, regardless of the advertising spend, the predicted revenue is always the same constant value.
import numpy as np
# Computing average revenue from the training set
average_revenue = y_train.mean()
# Calculating MAE using the baseline (mean prediction)
mae_baseline = np.mean(np.abs(y_test - average_revenue))
# Printing results
mae_baselinenp.float64(13580.480177514795)The MAE is now much higher. This shows that even a simple linear regression model can produce substantially better predictions than a naive baseline that ignores all features and uses only the average.
Multiple Linear Regression in Machine Learning
Up to now, we have not introduced any entirely new concept. Instead of creating a linear regression model using the whole original dataset, we split the dataset into two parts: one to estimate coefficients and the other to make predictions and assess accuracy. We continue the same approach here to see how including more predictors can affect prediction performance.
This time, we add the feature Product_Quality_Score to the model and train a multiple linear regression:
# Separating the predictors from the target variable
X = eshop_revenues[["Ad_Spend", "Product_Quality_Score"]]
y = eshop_revenues["Revenue"]
# Splitting data into training and test sets
X_train, X_test, y_train, y_test = train_test_split(
X, y, test_size = 0.20, random_state = 42
)
# Creating a linear regression model
lm_model_multiple = LinearRegression()
# Estimating the model parameters using the training data
lm_model_multiple.fit(X_train, y_train)We follow the same steps as before to calculate the MAE on the test set:
# Making predictions on the test_set
y_pred = lm_model_multiple.predict(X_test)
# Calculating the mean absolute error
mae = mean_absolute_error(y_test, y_pred)
# Printing the result
mae7472.2474603022865This time, the MAE is much lower! Adding the additional predictor improves the model’s predictions. Can we do even better? Let’s create a more advanced model by including both Ad_Spend and Product_Quality_Score along with their interaction:
# Creating interaction term manually
eshop_revenues["AdSpend_x_Quality"] = (
eshop_revenues["Ad_Spend"] * eshop_revenues["Product_Quality_Score"]
)
# Separating the predictors from the target variable
X = eshop_revenues[["Ad_Spend", "Product_Quality_Score", "AdSpend_x_Quality"]]
y = eshop_revenues["Revenue"]
# Splitting data into training and test sets
X_train, X_test, y_train, y_test = train_test_split(
X, y, test_size = 0.2, random_state = 42
)
# Creating a linear regression model
lm_model_multiple = LinearRegression()
# Estimating the model parameters using the training data
lm_model_multiple.fit(X_train, y_train)
# Making predictions on the test_set
y_pred = lm_model_multiple.predict(X_test)
# Calculating the mean absolute error
mae = mean_absolute_error(y_test, y_pred)
# Printing the result
mae7014.642195395383The MAE decreases even further, showing that including the interaction of the two predictors improves model performance.
Now, let’s push things to the extreme by including all available features and their pairwise interactions, and then evaluate the model’s accuracy on the test set. To do this, we first need to construct an expanded feature space that includes not only the original variables, but also all interaction terms between them. In scikit-learn, this can be achieved using the PolynomialFeatures class from the sklearn.preprocessing module, which automatically generates these additional features.
Before applying this transformation, we define X simply by dropping the target variable (Revenue) from the dataset along with the created variable AdSpend_x_Quality. This ensures that all remaining columns are treated as potential predictors, and that no information from the target variable leaks into the feature set:
from sklearn.preprocessing import PolynomialFeatures
# Defining predictors and target
X = eshop_revenues.drop(columns = ["Revenue", "AdSpend_x_Quality"])
y = eshop_revenues["Revenue"]We then apply PolynomialFeatures with degree = 2, which generates all pairwise interactions between features, in addition to the original variables. Since we are only interested in interaction effects and not higher-order polynomial terms, this step effectively expands our feature space to include all second-order relationships between predictors.
# Creating a polynomial feature generator to include all pairwise interactions
poly = PolynomialFeatures(degree = 2, interaction_only = True, include_bias = False)
# Expanding the feature matrix into polynomial feature space
X_poly = poly.fit_transform(X)The resulting matrix X_poly now contains the original features as well as all pairwise interaction terms between them. We can then proceed with the usual workflow by splitting the transformed dataset into training and test sets, fitting a linear regression model, and evaluating its predictive performance:
# Splitting data into training and test sets
X_train, X_test, y_train, y_test = train_test_split(
X_poly, y, test_size = 0.2, random_state = 42
)
# Creating a linear regression model
lm_model_all = LinearRegression()
# Estimating the model parameters using the training data
lm_model_all.fit(X_train, y_train)
# Making predictions on the test_set
y_pred = lm_model_all.predict(X_test)
# Calculating the mean absolute error
mae = mean_absolute_error(y_test, y_pred)
# Printing the result
mae5151.846199819311The performance looks solid. On average, predicted revenues deviate by roughly €5,150 from the actual revenues.
Does this mean we should always include all features and interactions? Not necessarily. The bias-variance trade-off reminds us that adding more features and interactions can sometimes increase variance and reduce generalization, leading to a worse test accuracy.
To see this in action, let’s follow the same process as before—use all available features and their interactions, but this time exclude the feature New_Product:
# Defining predictors and target
X = eshop_revenues.drop(columns = ["Revenue", "AdSpend_x_Quality", "New_Product"])
y = eshop_revenues["Revenue"]
# Creating a polynomial feature generator to include all pairwise interactions
poly = PolynomialFeatures(degree = 2, interaction_only = True, include_bias = False)
# Expanding the feature matrix into polynomial feature space
X_poly = poly.fit_transform(X)
# Splitting data into training and test sets
X_train, X_test, y_train, y_test = train_test_split(
X_poly, y, test_size = 0.2, random_state = 42
)
# Creating a linear regression model
lm_model_all_except_new = LinearRegression()
# Estimating the model parameters using the training data
lm_model_all_except_new.fit(X_train, y_train)
# Making predictions on the test_set
y_pred = lm_model_all_except_new.predict(X_test)
# Calculating the mean absolute error
mae = mean_absolute_error(y_test, y_pred)
# Printing the result
mae5020.158675621228Surprisingly, MAE is even lower—approximately €5,020—when we use a model with fewer predictors. This is a classic example of the bias-variance trade-off: a more complex model may fit the training data very well (low bias) but generalize worse (high variance).
We might think that this feature simply wasn’t important in the first place. Just like including irrelevant independent variables in a regression model for statistical inference can increase the standard errors of all coefficients without improving estimates (see Chapter Multiple Linear Regression and Statistical Inference), including such a feature in a machine learning model can hurt predictive performance.
But was New_Product really irrelevant? Let’s check the R-squared and the p-values of lm_model_all and lm_model_all_except_new. The code below initially uses the same machine learning workflow but then switches to a slightly different approach when finalizing the dataset with the predictor and target values. It then uses the statsmodels module to fit a linear regression model in order to extract the R-squared and the p-values of the two models. Our focus is not the code itself but on the statistical interpretation of the results and the comparison between the two model specifications.
import statsmodels.api as sm
from statsmodels.regression.linear_model import OLS
# Defining predictors and target
X = eshop_revenues.drop(columns=["Revenue", "AdSpend_x_Quality"])
y = eshop_revenues["Revenue"]
# Creating polynomial features (all pairwise interactions)
poly = PolynomialFeatures(degree = 2, interaction_only = True, include_bias = False)
X_poly = poly.fit_transform(X)
# Keep feature names (IMPORTANT for OLS interpretability)
feature_names = poly.get_feature_names_out(X.columns)
X_poly = pd.DataFrame(X_poly, columns=feature_names)
# Splitting data
X_train, X_test, y_train, y_test = train_test_split(
X_poly, y, test_size = 0.2, random_state = 42
)
# Reconstruct training dataset
train_data = X_train.copy()
train_data["Revenue"] = y_train
# Fitting linear regression model
lm_model_all_check = OLS(
train_data["Revenue"],
sm.add_constant(train_data.drop(columns = ["Revenue"]))
).fit()
# Printing p-values
lm_model_all_check.pvaluesconst 3.16e-02
Ad_Spend 1.96e-03
Website_Visitors 1.49e-01
Product_Quality_Score 7.95e-02
Customer_Satisfaction 1.58e-02
Price 5.64e-03
New_Product 3.77e-01
Ad_Spend Website_Visitors 1.51e-02
Ad_Spend Product_Quality_Score 2.08e-05
Ad_Spend Customer_Satisfaction 2.98e-02
Ad_Spend Price 3.52e-08
Ad_Spend New_Product 1.21e-01
Website_Visitors Product_Quality_Score 7.30e-01
Website_Visitors Customer_Satisfaction 1.00e-01
Website_Visitors Price 1.53e-01
Website_Visitors New_Product 5.89e-01
Product_Quality_Score Customer_Satisfaction 6.77e-03
Product_Quality_Score Price 2.14e-01
Product_Quality_Score New_Product 4.90e-02
Customer_Satisfaction Price 1.47e-02
Customer_Satisfaction New_Product 5.15e-02
Price New_Product 9.81e-01
dtype: float64# Printing R-squared
lm_model_all_check.rsquarednp.float64(0.9110229770490959)# Defining predictors and target
X = eshop_revenues.drop(columns=["Revenue", "AdSpend_x_Quality", "New_Product"])
y = eshop_revenues["Revenue"]
# Creating polynomial features (all pairwise interactions)
poly = PolynomialFeatures(degree = 2, interaction_only = True, include_bias = False)
X_poly = poly.fit_transform(X)
# Keep feature names (IMPORTANT for OLS interpretability)
feature_names = poly.get_feature_names_out(X.columns)
X_poly = pd.DataFrame(X_poly, columns = feature_names)
# Splitting data
X_train, X_test, y_train, y_test = train_test_split(
X_poly, y, test_size = 0.2, random_state = 42
)
# Reconstruct training dataset
train_data = X_train.copy()
train_data["Revenue"] = y_train
# Fitting linear regression model
lm_model_all_except_new_check = OLS(
train_data["Revenue"],
sm.add_constant(train_data.drop(columns = ["Revenue"]))
).fit()
# Printing R-squared
lm_model_all_except_new_check.rsquarednp.float64(0.8868067648123408)Looking at the summary of the model-fit, lm_model_all seemed actually better. Not only was R-Squared (goodness-of-fit metric) slightly higher, but there were also two marginally statistically significant interaction effects. However, the model lm_model_all_except_new had a better accuracy on the test set. This highlights an important lesson: a model optimized for statistical inference is not necessarily the same as one optimized for prediction. Different purposes may require different models.
Disclaimer
Creating and testing many different models and checking their accuracy on the test set is not best practice because it effectively makes the test set part of the training process. We discuss better approaches in the later chapters.
Recap
In this chapter, we explored simple and multiple linear regression from a predictive perspective, focusing on model evaluation using unseen data rather than statistical inference. We used the scikit-learn workflow to train models, generate predictions, and assess performance using MAE, gradually increasing model complexity by adding predictors, interaction terms, and full feature expansions. We showed that while more complex models often improve training fit, they do not always improve test performance, highlighting the bias–variance trade-off and the importance of evaluating models based on their predictive accuracy.