2.1. Linear Models

The simplest Machine Learning (ML) models are based on the assumption that the target can be expressed using a linear combination of the input features. That intuition leads to a class of models called linear models.

Linear Regression

Let \(\mathcal D = \{\boldsymbol x_i, y_i\}_{i=1}^N\) be a dataset with \(N\) observations and \(n\) features. The linear regression model describes the continuous target by a linear combination of the inputs, that is: \[ \boldsymbol{\hat{y}} = \boldsymbol{Ax}^T + b, \tag{1}\] where \(\boldsymbol{\hat{y}} \in \mathbb R^N\) is the vector with the predictions of the targets, \(\boldsymbol x \in \mathbb R^{N\times n}\) with the features of the observations, \(\boldsymbol A \in \mathbb R^n\) is the coefficient vector and \(b \in \mathbb R\) is the bias term. For each observation, the model leaves residuals of the form \(y_i - \hat{y}_i = \epsilon_i\). These are minimized using least-squares approaches, in order to find the best values for \(\boldsymbol A\) and \(b\).

Though simple, this model is intrinsically interpretable, as taking the partial derivatives with respect to each input feature yields: \[ \frac{\partial \boldsymbol{\hat y}}{\partial x_i} = a_i. \tag{2}\] The bigger the value \(a_i\), the more impact the feature \(x_i\) has on the outputs. Note that this could be influenced by the scale of the features. It is possible to add regularization to the linear regression, such as \(\ell_1\), known as Lasso regression, \(ell_2\), known as Ridge regression, or both, known as elastic net

Code: Linear Regression

Importing the libraries for a linear regression.

import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression

Generating the data.

X = np.random.uniform(-5, 5, 100).reshape(-1, 1)
y = 2.5 * X.squeeze() + 1.0 + np.random.normal(0, 2, 100)

Creating and fitting the model.

model = LinearRegression()
model.fit(X, y)

LinearRegression()

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

Printing the coefficients.

A = model.coef_[0]
b = model.intercept_

print(f"y_hat = {A:.2f} * x + {b:.2f}")

y_hat = 2.51 * x + 0.86

Plotting the final result.

plt.figure(figsize=(8, 5))
plt.scatter(X, y, color='blue', alpha=0.6, label='Observations')
plt.plot(X, model.predict(X), color='red', linewidth=2, label=f'y_hat = {A:.2f}x + {b:.2f}')
plt.xlabel("x")
plt.ylabel("y")
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)
plt.show()

Source: linear_regression.ipynb

Logistic Regression

The logistic regression is a classification model that measures the probability of a given data point belonging to a certain class. It utilizes a sigmoid function with a threshold value, typically \(P(\boldsymbol y_i=1|\boldsymbol x_i)>0.5\), for classification. The equation for logistic regression is given by: \[ P\left(\boldsymbol{y_i} = 1 \mid \boldsymbol{x}_i \right) = \frac{1}{1+e^{-\boldsymbol{A} \boldsymbol{x}_i^T}}. \tag{3}\]

Code: Logistic Regression

Importing the libraries for a logistic regression.

import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LogisticRegression

Generating the data.

np.random.seed(42)
X = np.random.normal(0, 2, 100).reshape(-1, 1)
y = (X[:, 0] > 0).astype(int)
flip_idx = np.random.choice(100, size=10, replace=False)
y[flip_idx] = 1 - y[flip_idx]

Creating and fitting the model.

model = LogisticRegression()
model.fit(X, y)

LogisticRegression()

Plotting the final result.

X_test = np.linspace(-6, 6, 300).reshape(-1, 1)
probabilities = model.predict_proba(X_test)[:, 1]

plt.figure(figsize=(8, 5))
plt.scatter(X, y, color='blue', alpha=0.5, label='Observations')
plt.plot(X_test, probabilities, color='red', linewidth=2, label=r'$P(y=1 | x)$')

plt.xlabel("x")
plt.ylabel("P")
plt.legend()
plt.show()

Source: logistic_regression.ipynb

Generalized Linear Models (GLMs)

Both linear and logistic regression are specific cases of a broader framework known as Generalized Linear Models (GLMs). GLMs assume that the target variable follows a distribution from the exponential family, and that its expected value is conditioned on the input features: \[ \mathbb{E}[y_i \mid \boldsymbol{x}_i] = \mu_i. \tag{4}\] To predict the targets, GLMs utilize a link function \(g(\cdot)\) that maps the expected value of the distribution to a linear predictor: \[ g(\mathbb{E}[y_i \mid \boldsymbol{x}_i]) = \boldsymbol{A}\boldsymbol{x}_i^T + b. \tag{5}\] By choosing specific families of distributions and link functions – such as an identity link for Gaussian targets or a logit link for binomial targets – the previous models are recovered.

Code: Generalized Linear Models

:::{.quarto-embed-nb-cell notebook=“/home/mzaiam/personal/docs/notes/machine_learning/supervised_methods/supervised_methods_assets/glm.ipynb” notebook-title=” Generalized Linear Model Regression Results ” notebook-cellId=“cell-0”}

Importing the libraries for a GLM.

import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt

Generating the data.

np.random.seed(42)
X = np.random.uniform(0, 3, 100)

mu = np.exp(1.2 * X + 0.5)
y = np.random.poisson(mu)

Creating and fitting the model.

X_sm = sm.add_constant(X)

glm_model = sm.GLM(y, X_sm, family=sm.families.Poisson(link=sm.families.links.Log()))
glm_results = glm_model.fit()

print(glm_results.summary())

                 Generalized Linear Model Regression Results                  
==============================================================================
Dep. Variable:                      y   No. Observations:                  100
Model:                            GLM   Df Residuals:                       98
Model Family:                 Poisson   Df Model:                            1
Link Function:                    Log   Scale:                          1.0000
Method:                          IRLS   Log-Likelihood:                -255.46
Date:                Fri, 24 Apr 2026   Deviance:                       120.93
Time:                        16:48:40   Pearson chi2:                     118.
No. Iterations:                     5   Pseudo R-squ. (CS):              1.000
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const          0.4624      0.088      5.240      0.000       0.289       0.635
x1             1.2107      0.038     31.939      0.000       1.136       1.285
==============================================================================

Plotting the final result.

X_test = np.linspace(0, 3, 100)
X_test_sm = sm.add_constant(X_test)
y_pred = glm_results.predict(X_test_sm)

plt.figure(figsize=(8, 5))
plt.scatter(X, y, color='blue', alpha=0.5, label='Observations')
plt.plot(X_test, y_pred, color='red', linewidth=2, label=r'GLM prediction')
plt.xlabel("x")
plt.ylabel("y")
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)
plt.show()

:::

Generalized Additive Models (GAMs)

Generalized Additive Models (GAMs) relax the strict linearity assumption of GLMs, allowing for more flexible, non-linear relationships while preserving interpretability. The general form of a GAM replaces the linear combination of features with a sum of smooth, feature-specific functions: \[ g(\mathbb{E}[y_i \mid \boldsymbol{x}_i]) = a_0 + \sum_{j=1}^n f_j(x_{i, j}), \tag{6}\] where \(f_j\) is an arbitrary function associated with the \(j\)-th feature. These functions are typically modeled using basis functions, such as splines. Because the model’s output is formulated as the sum of separate input functions, it remains highly interpretable. However, learning these separate non-linear functions without overfitting is a difficult task.

To address this, Explainable Boosting Machines (EBMs) (Lou et al. 2013) implement GAMs using cyclic gradient boosting. EBMs iteratively fit shallow decision trees to each feature independently, learning to correct the residual errors of previous trees. That is, at the \(i\)-th iteration, the tree associated to the \(j\)-th feature is trained to predict the residual \(r_{j-1}^{(i)}\), in a iterative way. As each tree is restricted to a single feature, it is considered a weak learner, that iteratively refines the corresponding function \(f_j\). After several iterations, the model produces an accurate yet transparent representation of the additive functions. Figure 1 depicts the diagram of an EBM.

Figure 1: Diagram of an EBM with \(n\) features and \(k\) iterations. In the end of the training, the trees of each feature are combined into a single linear function.

Code: Generalized Linear Models

Importing the libraries for an EBM.

import pandas as pd
import matplotlib.pyplot as plt
from sklearn.datasets import load_breast_cancer
from sklearn.model_selection import train_test_split
from interpret.glassbox import ExplainableBoostingClassifier

Loading a dataset.

data = load_breast_cancer()
X = pd.DataFrame(data.data, columns=data.feature_names)
y = data.target

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

Creating and fitting the model.

ebm = ExplainableBoostingClassifier(
    interactions=0,
    outer_bags=1,
    max_bins=64,
    max_rounds=1000,
    learning_rate=0.05,
    n_jobs=-1,
    random_state=42
)
ebm.fit(X_train, y_train)

ExplainableBoostingClassifier(interactions=0, learning_rate=0.05, max_bins=64,
                              max_rounds=1000, n_jobs=-1, outer_bags=1)

Plotting the contribution of the feature mean radius.

feature_name = 'mean radius'
feature_idx = X.columns.get_loc(feature_name)

ebm_global = ebm.explain_global()
feature_data = ebm_global.data(feature_idx)

x_vals = feature_data['names']
y_vals = feature_data['scores']

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

if len(x_vals) == len(y_vals) + 1:
    plt.step(x_vals[:-1], y_vals, where='post', color='blue', linewidth=2)
else:
    plt.plot(x_vals, y_vals, color='blue', linewidth=2)

if 'upper_bounds' in feature_data and 'lower_bounds' in feature_data:
    upper = feature_data['upper_bounds']
    lower = feature_data['lower_bounds']
    if len(x_vals) == len(upper) + 1:
        plt.fill_between(x_vals[:-1], lower, upper, alpha=0.2, color='blue', step='post')

plt.xlabel(feature_name)
plt.ylabel("Contribution")
plt.axhline(0, color='red', linestyle='--', alpha=0.5)
plt.grid(True, linestyle='--', alpha=0.7)
plt.show()

As depicted in the plot, lower values of the feature lead to higher probability of a benign tumor, while higher values increase the probability of malign tumors.

Source: ebm.ipynb

References

Lou, Yin, Rich Caruana, Johannes Gehrke, and Giles Hooker. 2013. “Accurate Intelligible Models with Pairwise Interactions.” In International Conference on Knowledge Discovery and Data Mining, 623–31.