Generalised Linear Models

Neil D. Lawrence

Back to lecture list

Generalised Linear Models

Week 6: Generalised Linear Models

[jupyter][google colab][reveal][edit]

Neil D. Lawrence

Abstract:

Setup

[edit]

notutils

[edit]

This small package is a helper package for various notebook utilities used below.

The software can be installed using

import importlib.util

from the command prompt where you can access your python installation.

The code is also available on GitHub: https://github.com/lawrennd/notutils

Once notutils is installed, it can be imported in the usual manner.

import notutils

pods

[edit]

In Sheffield we created a suite of software tools for ‘Open Data Science.’ Open data science is an approach to sharing code, models and data that should make it easier for companies, health professionals and scientists to gain access to data science techniques.

You can also check this blog post on Open Data Science.

The software can be installed using

cmd = install_command('pods')

%system {cmd}

from the command prompt where you can access your python installation.

The code is also available on GitHub: https://github.com/lawrennd/ods

Once pods is installed, it can be imported in the usual manner.

import pods

mlai

[edit]

The mlai software is a suite of helper functions for teaching and demonstrating machine learning algorithms. It was first used in the Machine Learning and Adaptive Intelligence course in Sheffield in 2013.

The software can be installed using

cmd = install_command('mlai')

%system {cmd}

from the command prompt where you can access your python installation.

The code is also available on GitHub: https://github.com/lawrennd/mlai

Once mlai is installed, it can be imported in the usual manner.

import mlai
from mlai import plot

cmd = install_command('statsmodels')

%system {cmd}

Review

We introduced machine learning as a way to extract knowledge from data to make predictions through a prediction function and an objective function. We looked at a simple example of predicting whether someone would buy a jumper based on their age and latitude, using logistic regression to model the log-odds of purchase. This highlighted how machine learning can codify predictions through mathematical functions. This is an example of a broader approach known as generalized linear models.

When taking a probabilistic approach to supervised learning we’re interested in predicting a class label, $y_i$, given an input, $\mathbf{ x}_i$. That’s represented probabilisticially as $p(y_i|\mathbf{ x}_i)$. We can derive this conditional distribution through either (1) modelling the joint distribution, $p(\mathbf{ y}, \mathbf{X})$ and then dividing by the marginal distribution of the inputs, $p(\mathbf{X})$ , or (2) focusing specifically on modeling the conditional density, $p(\mathbf{ y}|\mathbf{X})$, that directly answers our prediction question. In the generalised linear model we choose the second approach.

As we move to generalized linear models like logistic regression, we’ll see how directly modeling the conditional density $p(\mathbf{ y}|\mathbf{X})$ can provide more flexibility in our modeling assumptions, while still allowing us to make the specific predictions we need.

Model Diagnostics and Linear Regression

[edit]

In linear regression, we model the relationship between a continuous response variable $y_i$ and input variables $\mathbf{ x}_i$ through a linear function with Gaussian noise:

\[y_i = f(\mathbf{ x}_i) + \epsilon_i\]

where $f(\mathbf{ x}_i) = \mathbf{ w}^\top\mathbf{ x}_i = \sum_{j=1}^D w_jx_{i,j}$ and $\epsilon_i \sim \mathscr{N}\left(0,\sigma^2\right)$

This gives us a probabilistic model:

\[p(y_i|\mathbf{ x}_i) = \mathscr{N}\left(y_i|\mathbf{ w}^\top\mathbf{ x}_i,\sigma^2\right)\]

The key components are:

$y_i$ is the target/response variable we want to predict
$\mathbf{ x}_i$ contains the input features/explanatory variables
$\mathbf{ w}$ contains the parameters/coefficients we learn
$\epsilon_i$ represents random Gaussian noise with variance $\sigma^2$

For the full dataset, we can write this in matrix form:

\[\mathbf{ y}= \mathbf{X}\mathbf{ w}+ \boldsymbol{ \epsilon}\]

where $\mathbf{ y}= [y_1,\ldots,y_N]^\top$, $\mathbf{X}$ contains the input vectors as rows, and $\boldsymbol{ \epsilon}\sim \mathscr{N}\left(\mathbf{0},\sigma^2\mathbf{I}\right)$.

The expected value of our prediction is:

\[\mathbb{E}[y_i|\mathbf{ x}_i] = \mathbf{ w}^\top\mathbf{ x}_i\]

This linear model forms the foundation for generalized linear models like logistic regression, where we’ll adapt the model for classification by transforming the output through a non-linear function.

import statsmodels.api as sm
import pods

# Demo of linear regression using python statsmodels.
data = pods.datasets.olympic_marathon_men()
xlim = (1876,2044)
x = data['X']
y = data['Y']

The statsmodels package makes it easy to add the constant term to the matrix.

# Add constant term to design matrix
print(x.shape)
x = sm.add_constant(x)
print(x.shape)

model = sm.OLS(y, x)

Using statsmodels we don’t need to implement the QR decomposition directly we can simply fit the model.

results = model.fit()

The model results also contain a lot more than just the regression weights. Various diagnostic statistics are also computed including checks on the residuals (i.e. $y_i - f_i$) which should be Gaussian (or normal) distributed.

We can show obtain our predictions, $f$ from the model as follows.

f = results.predict(x)

And plot them against the actual data.

Figure: Plot of the statsmodels fit to the olympic marathon data.

results.summary()

The summary helps us evaluate our model fit and identify potential areas for improvement. Since we’re working with one-dimensional data (year vs time), we can visualise everything easily to complement these statistical measures.

The model fit statistics show a moderately strong fit, with an R-squared of 0.730 indicating that our model explains 73.0% of the variance in the data. The adjusted R-squared of 0.721 confirms this isn’t just due to overfitting (which is unsurprising given we are only fitting a linear model). The very low F-statistic p-value (1.85e-09) confirms the model’s overall significance.

The AIC (11.33) and BIC (14.13) information criteria values are used to compare the model against alternatives we might try.

Looking at the model parameters, we see a coefficient of -0.0116 for our predictor, with a small standard error (0.001). The $t$-statistic of -8.710 and p-value of less than 0.001 indicate this effect is highly significant. The 95% confidence interval [-0.014, -0.009] gives us good confidence in our estimate. The negative coefficient confirms the expected downward trend in marathon times over the years.

Figure: The histogram of the residuals of the fit to the Olympic Marathon data.

We can see in the histogram that the residuals are not Gaussian (normal) distributed. Looking at the diagnostics we see the same challenges. The histogram doesn’t show us the autocorrelation, that comes from the Durbin-Watson statistic.

The Durbin-Watson statistic (0.981) indicates positive autocorrelation in the residuals. This suggests we might want to
- Consider time series modeling approaches
- Add polynomial terms to capture non-linear trends
- Investigate if there are distinct “eras” in marathon times
The Jarque-Bera test ($p$-value 8.32e-14) tells us our residuals aren’t normally distributed. The skew (1.947) and kurtosis (8.746) values show the distribution is strongly right-skewed with very heavy tails. This is due to the outlier in 1904. We could deal with this through transformation of the response variable, robust regression or by dropping the outlier.
The large condition number (9.94e+04) suggests potential numerical instability or multicollinearity issues. This is because we haven’t centred and scaled the model.

The beauty of having one-dimensional data is that we can plot everything to visually confirm these statistical findings. It allows us to visually assess the linearity assumption, identify potential outliers (such as the 1904 point) and spot any systematic patterns in the residuals.

This visual inspection, combined with our statistical diagnostics, will guide our next steps in improving the model.

Figure: Linear regression fit to Olympic marathon men’s times using statsmodels.

Looking at our plot and model diagnostics, we can now better understand the large condition number (1.08e+05) in our results. This high value likely stems from using raw year values (e.g., 1896, 1900, etc.) as our predictor variable. Such large numbers can lead to numerical instability in the computations.

To address this, we will look at centering the years around their mean and scaling the years to a smaller range.

Looking at the plot we can see more about our our data. The 1904 St. Louis Olympics appears as a clear outlier, contributing to the non-normal residuals (Jarque-Bera p=0.00432) and right-skewed distribution (skew=1.385).

There are also regimes in the data, such as the rapid improvement in times pre-WWI, the disruption and variation during the war years, the more steady, consistent progress post-WWII. These regime changes help explain the strong positive autocorrelation (Durbin-Watson=0.242) in our residuals.

This suggests we could improve the model by adding additional features such as polynomial terms to capture non-linear trends, indicator variables for different time periods and interaction terms between features. We could also consider variables accounting for external factors like temperature or course conditions.

To incorporate multiple features into our model, we need a systematic way to organize this additional information. This brings us back to the concept of the design matrix.

Design Matrix

The design matrix, often denoted as $\boldsymbol{ \Phi}$, is a key component of a statistical model. It organizes our feature data in a structured way that facilitates model fitting and analysis. Each row of the design matrix represents a single observation or data point, while each column represents a different feature or predictor variable.

For $n$ observations and $p$ features, the design matrix takes the form:

\[\boldsymbol{ \Phi}= \begin{bmatrix} \phi_{11} & \phi_{12} & \cdots & \phi_{1p} \\ \phi_{21} & \phi_{22} & \cdots & \phi_{2p} \\ \vdots & \vdots & \ddots & \vdots \\ \phi_{n1} & \phi_{n2} & \cdots & \phi_{np} \end{bmatrix}\]

For example, if we’re predicting house prices, each row might represent a different house, with columns for features like: square footage, number of bedrooms, yyear built, plot size.

The design matrix provides a compact way to represent all our feature data and is used directly in model fitting. When we write our linear model as $\mathbf{ y}= \boldsymbol{ \Phi}\mathbf{ w}+ \boldsymbol{ \epsilon}$, the design matrix $\boldsymbol{ \Phi}$ is multiplied by our parameter vector $\mathbf{ w}$ to generate predictions.

import statsmodels.api as sm
import pods
import numpy as np

# Demo of additional features with interactions regression usying python statsmodels.
data = pods.datasets.olympic_marathon_men()
x = data['X']
y = data['Y']
xlim = (1876,2044)

Now we process our input data to create the design matrix.

def marathon_design_matrix(x, product_terms=True):
    """This function augments the input matrix with indicator variables that state the year."""
    
    # Scale the year to avoid numerical issues
    x_scaled = (x - 1900) / 100  # Center around 1900 and scale to century units
    Phi = x_scaled

    # Add to design matrix indicator variable for 1914-1945
    Phi = np.hstack([Phi, ((x[:, 0] >= 1914) & (x[:, 0] <= 1945)).astype(np.float64)[:, np.newaxis]])

    # Add to design matrix indicator variable for post-1945
    Phi = np.hstack([Phi, (x[:, 0] > 1945).astype(np.float64)[:, np.newaxis]])
    if product_terms:
        # Add product terms that multiply the scaled year and the indicator variables.
        Phi = np.hstack([Phi, x_scaled[:, 0:1] * Phi[:, 1:2], x_scaled[:, 0:1] * Phi[:, 2:3]])
    Phi = sm.add_constant(Phi)
    return Phi

# Create the model
Phi = marathon_design_matrix(x)
model = sm.OLS(y, Phi)

Now we can fit the model.

# Do the linear fit
results = model.fit()

Figure: Plot of the fit obtained when we include indicator variables that identify the pre-1914, across war, and post-war eras.

Let’s check the summary of the results.

results.summary()

Using indicator variables for historical periods (and without relying on quadratic terms), the model shows a strong improvement over the simpler linear model, with an R-squared of 0.877 and an adjusted R-squared of 0.851. This indicates that about 87% of the variance in marathon times is explained by our model. The overall model significance is high (F-statistic $p=3.70\times 10^{-10}$), and the information criteria are AIC $=-4.164$, BIC $=4.244$.

Key parameter estimates (with 95% confidence intervals) are:

Intercept (const): $4.5876$ (SE $0.113$), $t=40.533$, $p=0.000$, CI $[4.354, 4.821]$
$x_1$ (base time trend): $-3.6127$ (SE $1.634$), $t=-2.211$, $p=0.037$, CI $[-6.984, -0.241]$
$x_2$ (war period indicator): $-0.7523$ (SE $0.480$), $t=-1.567$, $p=0.130$, CI $[-1.743, 0.239]$
$x_3$ (post war indicator): $-0.9399$ (SE $0.211$), $t=-4.449$, $p=0.000$, CI $[-1.376, -0.504]$
$x_4$ (wartime time trend): $2.8799$ (SE $2.310$), $t=1.247$, $p=0.225$, CI $[-1.888, 7.648]$
$x_5$ (post war time trend): $3.0493$ (SE $1.646$), $t=1.853$, $p=0.076$, CI $[-0.348, 6.446]$

These results suggest a strong negative time trend and statistically significant effects for some period indicators, with others being suggestive but not conclusive.

However, there are some concerns:

The Jarque–Bera test ($p=1.34\times 10^{-24}$) indicates non-normal residuals. Skewness is $2.393$ and kurtosis is $11.065$.
The Durbin–Watson statistic is $2.101$, suggesting little residual autocorrelation.
The condition number is $127$, indicating the design is reasonably well conditioned relative to polynomial/interaction-augmented alternatives.

The indicator-variable model captures key regime differences and substantially improves fit over the baseline linear model. One side-effect of using abrupt indicator boundaries is that the fitted curve is discontinuous at regime change years. This can create edge effects where adjacent observations on either side of a boundary receive different offsets even if they are close in time.

If continuity is needed we could consider

Replacing hard indicators with smooth transition functions (e.g., logistic ramps centered on boundary years)
Adding interaction terms with polynomial or spline bases to allow the trend to bend smoothly through transitions
Using piecewise-linear or spline regression with continuity (and optionally slope) constraints at knots
Hierarchical/regularized models that shrink jump sizes toward zero unless strongly supported by data

These approaches preserve the idea of regime structure while avoiding artificial jumps at the boundaries.

Logistic Regression

[edit]

A logistic regression is an approach to classification which extends the linear regression models we’ve already explored. Rather than modeling the output of the function directly the assumption is that we model the log-odds with the basis functions.

The odds are defined as the ratio of the probability of a positive outcome, to the probability of a negative outcome. Just as we saw in our jumper (sweater) example where \[ \log \frac{p(\text{bought})}{p(\text{not bought})} = w_0 + w_1 \text{age} + w_2 \text{latitude} \] If the probability of a positive outcome is denoted by $\pi$, then the odds are computed as $\frac{\pi}{1-\pi}$.

Odds are widely used by bookmakers in gambling, although a bookmakers odds won’t normalise: i.e. if you look at the equivalent probabilities, and sum over the probability of all outcomes the bookmakers are considering, then you won’t get one. This is how the bookmaker makes a profit. Because a probability is always between zero and one, the odds are always between $0$ and $\infty$. If the positive outcome is unlikely the odds are close to zero, if it is very likely then the odds become close to infinite. Taking the logarithm of the odds maps the odds from the positive half space to being across the entire real line. Odds that were between 0 and 1 (where the negative outcome was more likely) are mapped to the range between $-\infty$ and $0$. Odds that are greater than 1 are mapped to the range between $0$ and $\infty$. Considering the log odds therefore takes a number between 0 and 1 (the probability of positive outcome) and maps it to the entire real line. The function that does this is known as the logit function, $g(p_i) = \log\frac{p_i}{1-p_i}$. This function is known as a link function.

For a standard regression we take, \[ f(\mathbf{ x}) = \mathbf{ w}^\top \boldsymbol{ \phi}(\mathbf{ x}), \] if we want to perform classification we perform a logistic regression. \[ \log \frac{\pi}{(1-\pi)} = \mathbf{ w}^\top \boldsymbol{ \phi}(\mathbf{ x}) \] where the odds ratio between the positive class and the negative class is given by \[ \frac{\pi}{(1-\pi)} \] The odds can never be negative, but can take any value from 0 to $\infty$. We have defined the link function as taking the form $g^{-1}(\cdot)$ implying that the inverse link function is given by $g(\cdot)$. Since we have defined, \[ g(\pi) = \mathbf{ w}^\top \boldsymbol{ \phi}(\mathbf{ x}) \] we can write $\pi$ in terms of the inverse link function, $h(\cdot)$ as \[ \pi = h(\mathbf{ w}^\top \boldsymbol{ \phi}(\mathbf{ x})). \]

We’ll define our prediction, objective and gradient functions below. But before we start, we need to define a basis function for our model. Let’s start with the linear basis.

import numpy as np

from mlai import linear

Soft Arg Max (Softmax) Function

[edit]

Given class scores $f_k(\mathbf{ x})=\mathbf{ w}_k^\top\boldsymbol{ \phi}(\mathbf{ x})$, the softmax maps scores to probabilities: $\pi_k = \frac{\exp(f_k)}{\sum_{j=1}^K \exp(f_j)}$.

Prediction Function

[edit]

Now we have the basis function let’s define the prediction function.

import numpy as np

def predict(w, x, basis=mlai.linear, **kwargs):
    "Generates the prediction function and the basis matrix."
    Phi = basis(x, **kwargs)
    f = np.dot(Phi, w)
    return 1./(1+np.exp(-f)), Phi

This inverse of the link function is known as the logistic (thus the name logistic regression) or sometimes it is called the sigmoid function. For a particular value of the input to the link function, $f_i = \mathbf{ w}^\top \boldsymbol{ \phi}(\mathbf{ x}_i)$ we can plot the value of the inverse link function as below.

By replacing the inverse link with the sigmoid we can write $\pi$ as a function of the input and the parameter vector as, \[ \pi(\mathbf{ x},\mathbf{ w}) = \frac{1}{1+\exp\left(-\mathbf{ w}^\top \boldsymbol{ \phi}(\mathbf{ x})\right)}. \] The process for logistic regression is as follows. Compute the output of a standard linear basis function composition ($\mathbf{ w}^\top \boldsymbol{ \phi}(\mathbf{ x})$, as we did for linear regression) and then apply the inverse link function, $g(\mathbf{ w}^\top \boldsymbol{ \phi}(\mathbf{ x}))$. In logistic regression this involves squashing it with the logistic (or sigmoid) function. Use this value, which now has an interpretation as a probability in a Bernoulli distribution to form the likelihood. Then we can assume conditional independence of each data point given the parameters and develop a likelihod for the entire data set.

The Bernoulli likelihood is of the form, \[ P(y_i|\mathbf{ w}, \mathbf{ x}) = \pi_i^{y_i} (1-\pi_i)^{1-y_i} \] which we can think of as clever trick for mathematically switching between two probabilities if we were to write it as code it would be better described as

def bernoulli(x, y, pi):
    if y == 1:
        return pi(x)
    else:
        return 1-pi(x)

but writing it mathematically makes it easier to write our objective function within a single mathematical equation.

Maximum Likelihood

[edit]

To obtain the parameters of the model, we need to maximise the likelihood, or minimise the objective function, normally taken to be the negative log likelihood. With a data conditional independence assumption the likelihood has the form, \[ P(\mathbf{ y}|\mathbf{ w}, \mathbf{X}) = \prod_{i=1}^nP(y_i|\mathbf{ w}, \mathbf{ x}_i). \] which can be written as a log likelihood in the form \[ \log P(\mathbf{ y}|\mathbf{ w}, \mathbf{X}) = \sum_{i=1}^n\log P(y_i|\mathbf{ w}, \mathbf{ x}_i) = \sum_{i=1}^n y_i \log \pi_i + \sum_{i=1}^n(1-y_i)\log (1-\pi_i) \] and if we take the probability of positive outcome for the $i$th data point to be given by \[ \pi_i = h\left(\mathbf{ w}^\top \boldsymbol{ \phi}(\mathbf{ x}_i)\right), \] where $h(\cdot)$ is the inverse link function, then this leads to an objective function of the form, \[ E(\mathbf{ w}) = - \sum_{i=1}^ny_i \log h\left(\mathbf{ w}^\top \boldsymbol{ \phi}(\mathbf{ x}_i)\right) - \sum_{i=1}^n(1-y_i)\log \left(1-h\left(\mathbf{ w}^\top \boldsymbol{ \phi}(\mathbf{ x}_i)\right)\right). \]

Negative Log Likelihood

[edit]

import numpy as np

def objective(h, y):
    "Computes the objective function."
    labs = np.asarray(y, dtype=float).flatten()
    posind = np.where(labs==1)
    negind = np.where(labs==0)
    return -np.log(h[posind, :]).sum() - np.log(1-h[negind, :]).sum()

As normal, we would like to minimise this objective. This can be done by differentiating with respect to the parameters of our prediction function, $\pi(\mathbf{ x};\mathbf{ w})$, for optimisation. The gradient of the likelihood with respect to $\pi(\mathbf{ x};\mathbf{ w})$ is of the form, \[ \frac{\text{d}E(\mathbf{ w})}{\text{d}\mathbf{ w}} = -\sum_{i=1}^n \frac{y_i}{h\left(\mathbf{ w}^\top \boldsymbol{ \phi}(\mathbf{ x})\right)}\frac{\text{d}h(f_i)}{\text{d}f_i} \boldsymbol{ \phi}(\mathbf{ x}_i) + \sum_{i=1}^n \frac{1-y_i}{1-h\left(\mathbf{ w}^\top \boldsymbol{ \phi}(\mathbf{ x})\right)}\frac{\text{d}h(f_i)}{\text{d}f_i} \boldsymbol{ \phi}(\mathbf{ x}_i) \] where we used the chain rule to develop the derivative in terms of $\frac{\text{d}h(f_i)}{\text{d}f_i}$, which is the gradient of the inverse link function (in our case the gradient of the sigmoid function).

So the objective function now depends on the gradient of the inverse link function, as well as the likelihood depends on the gradient of the inverse link function, as well as the gradient of the log likelihood, and naturally the gradient of the argument of the inverse link function with respect to the parameters, which is simply $\boldsymbol{ \phi}(\mathbf{ x}_i)$.

The only missing term is the gradient of the inverse link function. For the sigmoid squashing function we have, \[\begin{aligned} h(f_i) &= \frac{1}{1+\exp(-f_i)}\\ &=(1+\exp(-f_i))^{-1} \end{align*}$$ and the gradient can be computed as $$\begin{align*} \frac{\text{d}h(f_i)}{\text{d} f_i} & = \exp(-f_i)(1+\exp(-f_i))^{-2}\\ & = \frac{1}{1+\exp(-f_i)} \frac{\exp(-f_i)}{1+\exp(-f_i)} \\ & = h(f_i) (1-g(f_i)) \end{aligned}\]

so the full gradient can be written down as \[ \frac{\text{d}E(\mathbf{ w})}{\text{d}\mathbf{ w}} = -\sum_{i=1}^n y_i\left(1-h\left(\mathbf{ w}^\top \boldsymbol{ \phi}(\mathbf{ x})\right)\right) \boldsymbol{ \phi}(\mathbf{ x}_i) + \sum_{i=1}^n (1-y_i)\left(h\left(\mathbf{ w}^\top \boldsymbol{ \phi}(\mathbf{ x})\right)\right) \boldsymbol{ \phi}(\mathbf{ x}_i). \]

import numpy as np

def gradient(h, Phi, y):
    "Generates the gradient of the parameter vector."
    labs = np.asarray(y, dtype=float).flatten()
    posind = np.where(labs==1)
    dw = -(Phi[posind]*(1-h[posind])).sum(0)
    negind = np.where(labs==0 )
    dw += (Phi[negind]*h[negind]).sum(0)
    return dw[:, None]

Optimisation of the Function

Reorganising the gradient to find a stationary point of the function with respect to the parameters $\mathbf{ w}$ turns out to be impossible. Optimisation has to proceed by numerical methods. Options include the multidimensional variant of Newton’s method or gradient based optimisation methods like we used for optimising matrix factorisation for the movie recommender system. We recall from matrix factorisation that, for large data, stochastic gradient descent or the Robbins Munro (Robbins and Monro, 1951) optimisation procedure worked best for function minimisation.

Multiplying everything out leads to \[ \begin{aligned} \frac{\text{d}E(\mathbf{ w})}{\text{d}\mathbf{ w}} = & -\sum_{i=1}^n y_i- h\left(\mathbf{ w}^\top \boldsymbol{ \phi}(\mathbf{ x})\boldsymbol{ \phi}(\mathbf{ x}_i). \]

Ad Matching for Facebook

[edit]

Logistic regression is a widely used technique in many real world application domains. Within social media, the first wave of machine learning solutions were targeted at challenges such as ad matching or news-post ranking. In these domains the users are characterised by their previous behaviour (such as which posts they have liked) and other contextual information (such as which friends they have). The logistic regression approach is then used to predict which ads should be shown. For example Facebook’s ad matching system faces the complex challenge of connecting millions of advertisers with billions of users. The scale of this matching problem is immense:

Facebook has millions of potential advertisers, each with different target audiences and budgets
There are billions of users, each with unique interests and behaviors
Decisions about which ads to show must be made in real-time
The system needs to balance between user experience and advertiser ROI

Facebook’s solution employed a system that combines logistic regression with decision trees. The system considers various features including:

User demographics and behaviour
Historical ad performance
Advertiser parameters and budgets
Contextual information

The hybrid approach allows Facebook to make rapid predictions about ad performance while handling the massive scale of their platform. The details of this system are described in a technical paper by Herbrich et al., which describes Facebook’s click prediction system for online advertising, which combines decision trees with logistic regression to predict ad clicks. At the time Facebook had over 750 million daily active users and 1 million advertisers at the time. The authors found that a hybrid model outperformed either decision trees or logistic regression alone.

Importantly, from our perspective, the researchers found that feature selection was the most critical factor - particularly choosing the features that capture historical information about users and ads. Given the ‘right’ features and model architecture, adapting the system to accommodate other aspects such as data freshness and learning rate tuning provided smaller incremental gains. To handle the scale of the task the system used a cascade of classifies to handle thecandidate ads per request efficiently. The hybrid model served as the final stage classifier. The paper also discusses practical considerations around implementing online learning and managing computational resources at scale.

Olivetti Glasses Data

[edit]

Let’s classify images with logistic regression. We’ll look at a data set of individuals with glasses. We can load in the data from pods as

Correspond to whether the subject of the image is wearing glasses. Set up the ipython environment and download the data:

from scipy import io

First let’s retrieve some data. We will use the ORL faces data set, our objective will be to classify whether or not the subject in an image is wearing glasess.

Here’s a simple way to visualise the data. Each pixel in the image will become an input to the GP.

import pods

data = pods.datasets.olivetti_glasses()
X = data['X']
y = data['Y']
Xtest = data['Xtest']
ytest = data['Ytest']
print(data['info'], data['details'], data['citation'])

Figure: Image from the Oivetti glasses data sets.

Batch Gradient Descent

[edit]

We will need to define some initial random values for our vector and then minimize the objective by descending the gradient.

# Separate train and test
indices = np.random.permutation(X.shape[0])
num_train = int(np.ceil(X.shape[0]/2))
train_indices = indices[:num_train]
test_indices = indices[num_train:]
X_train = X.iloc[train_indices]
y_train = y.iloc[train_indices]==True
X_test = X.iloc[test_indices]
y_test = y.iloc[test_indices]==True

import numpy as np

# gradient descent algorithm
w_vals = np.random.normal(size=(X.shape[1]+1, 1), scale = 0.001)
# Convert weights to a pandas DataFrame with column names
if typeof(X_train) is pd.DataFrame:
    w = pd.DataFrame(data=w_vals, 
                     index=['Eins'] + list(X_train.columns), 
                     columns=['weight'])

eta = 1e-9
iters = 10000
for i in range(iters):
    g, Phi = predict(w, X_train, mlai.linear)
    w -= eta*gradient(g, Phi, y_train) + 0.001*w
    if not i % 100:
        print("Iter", i, "Objective", objective(g, y_train))

Going Further: Optimization

[edit]

Other optimization techniques for generalized linear models include Newton’s method, it requires you to compute the Hessian, or second derivative of the objective function.

Methods that are based on gradients only include L-BFGS and conjugate gradients. Can you find these in python? Are they suitable for very large data sets? }

Other GLMs

[edit]

We’ve introduced the formalism for generalised linear models. Have a think about how you might model count data using the Poisson distribution and a log link function for the rate, $\lambda(\mathbf{ x})$. If you want a data set you can try the pods.datasets.google_trends() for some count data.

Poisson Distribution

[edit]

Figure: The Poisson distribution.

Poisson Regression

[edit]

Poisson regression is a type of generalized linear model (GLM) used when modeling count data. It assumes the response variable follows a Poisson distribution and uses a logarithmic link function to relate the mean of the response to the linear predictor.

In this model, we make the rate parameter λ a function of covariates (like space or time). The logarithm of the rate is modeled as a linear combination of the input features:

\[\log \lambda(\mathbf{ x}, t) = \mathbf{ w}_x^\top \boldsymbol{ \phi}_x(\mathbf{ x}) + \mathbf{ w}_t^\top \boldsymbol{ \phi}_t(t)\]

where: - $\mathbf{ w}_x$ and $\mathbf{ w}_t$ are parameter vectors - $\boldsymbol{ \phi}_x(\mathbf{ x})$ and $\boldsymbol{ \phi}_t(t)$ are basis functions for space and time respectively

This formulation is known as a log-linear or log-additive model because we’re adding terms in the log space. The logarithm serves as our link function, connecting the linear predictor to the response variable’s mean.

An important characteristic of this model that practitioners should be aware of is that while we add terms in the log space, the model becomes multiplicative when we transform back to the original space. This happens because:

We start with the log-additive form: $\log \lambda(\mathbf{ x}, t) = f_x(\mathbf{ x}) + f_t(t)$
When we exponentiate both sides to get back to λ, the addition in log space becomes multiplication: \[\lambda(\mathbf{ x}, t) = \exp(f_x(\mathbf{ x}) + f_t(t)) = \exp(f_x(\mathbf{ x}))\exp(f_t(t))\]

This multiplicative nature has important implications for interpretation. For example, if we increase one input variable, it has a multiplicative effect on the rate, not an additive one. This can lead to rapid growth in the predicted counts as input values increase.

Let’s look at another example using synthetic data to demonstrate Poisson regression without relying on external APIs.

import numpy as np
import statsmodels.api as sm

# Generate some example count data
np.random.seed(42)
n_samples = 100
x1 = np.random.uniform(0, 10, n_samples)
x2 = np.random.uniform(0, 5, n_samples)
X = np.column_stack((x1, x2))

# True relationship: y ~ Poisson(exp(1 + 0.3*x1 - 0.2*x2))
lambda_true = np.exp(1 + 0.3*x1 - 0.2*x2)
y = np.random.poisson(lambda_true)

# Fit Poisson regression
model_synthetic = sm.GLM(y, sm.add_constant(X), family=sm.families.Poisson())
result_synthetic = model_synthetic.fit()

The statsmodels library in Python provides a convenient way to fit Poisson regression models. The sm.GLM function is used to fit generalized linear models, and we specify the Poisson family to indicate that we’re modeling count data. The sm.add_constant(x) function adds a column of ones to the design matrix to account for the intercept term.

In this synthetic example, we generate count data that follows a Poisson distribution where the rate parameter $\lambda$ depends on two predictor variables. This demonstrates how Poisson regression can model count data with multiple predictors.

Figure: Diagnostic plots for the Poisson regression model showing actual vs predicted counts and residual analysis.

Practical Tips

[edit]

{In practice feature engineering is critical to the success of a GLM. Build modular data processing pipelines that allow you to easily test different feature sets. For example, if modeling house prices, you might want to test combinations of raw features (square footage, bedrooms), derived features (price per square foot), and interaction terms (bedrooms $\times$ bathrooms). Consider non-linear transformations of continuous variables. For instance, taking the log of price data often helps normalise distributions. Be thoughtful about encoding categorical variables. One-hot encoding isn’t always optimal, for high-cardinality categories, consider target encoding or feature hashing. Scale features appropriately: standardization or min-max scaling depending on your model assumptions. Document your feature creation process thoroughly, including the rationale for each transformation. Jupyter notebooks can be a good place to do this.

Model validation requires careful consideration. Cross-validation should match your real-world use case. For time series data, use time-based splits rather than random splits. Bootstrap sampling can help with understanding parameter uncertainty. For example, bootstrapping can show if a coefficient’s sign might flip under different samples. Hold-out test sets should be truly independent. In a customer churn model, this might mean testing on future customers rather than a random subset. Watch for data leakage, especially with time-dependent features. E.g. if predicting customer churn, using future purchase data would create leakage.

Diagnostic checks are essential for building confidence about the model reliability. Create residual plots against fitted values and each predictor. Look for systematic patterns, a U-shaped residual plot suggests missing quadratic terms. For logistic regression, plot predicted probabilities against actual outcomes in bins to check calibration. Calculate influence measures like Cook’s distance to identify outliers. In a house price model, a mansion might have outsized influence on coefficients. Check variance inflation factors (VIF) for multicollinearity. High VIF (>5-10) suggests problematic correlation between predictors.

Visualisation remains crucial throughout. Before modeling, create scatter plots, box plots, and histograms to understand your data distribution and relationships. Use pairs plots to identify correlations and potential interactions between features. Create residual diagnostic plots including Q-Q plots for normality checking. When communicating results, focus on interpretable visualizations. For instance, partial dependence plots can show how predictions change with a single feature.

Other practical considerations include starting simple and add complexity incrementally. A basic linear model often provides a good baseline. Keep track of model performance metrics across iterations to ensure changes actually improve results. Consider the computational cost of feature engineering: some transformations might not be feasible in production. Think about how features will be available in production. If a feature requires complex processing or external data, it might not be practical. For categorical variables with many levels, consider grouping rare categories. When dealing with missing data, document your imputation strategy and test its impact on model performance. When using the Access, Assess, Address approach remember that documentation of the missing values is part of Assess. Imputation can only occur in Address.

Thanks!

For more information on these subjects and more you might want to check the following resources.

References

Robbins, H., Monro, S., 1951. A stochastic approximation method. Annals of Mathematical Statistics 22, 400–407.

Rogers, S., Girolami, M., 2011. A first course in machine learning. CRC Press.

Generalised Linear Models

Week 6: Generalised Linear Models

Setup

notutils

pods

mlai

Review

Model Diagnostics and Linear Regression

Design Matrix

Logistic Regression

Soft Arg Max (Softmax) Function

Prediction Function

Maximum Likelihood

Negative Log Likelihood

Optimisation of the Function

Ad Matching for Facebook

Olivetti Glasses Data

Batch Gradient Descent

Going Further: Optimization

Other GLMs

Poisson Distribution

Poisson Regression

Practical Tips

Further Reading

Thanks!

References