Logistic regression

3 minute read

Published:

In this short tutorial we will implement two different methods for finding the maximum likelihood estimate (MLE) of a logistic regression model (McCullagh and Nelder, 1989), using Python. The MLE will be obtained using the algorithms:

  • Iteratively Reweighted Least Squares (IRLS)
  • Expectation Maximization (EM) via Polya-Gamma data augmentation

We start by loading some useful Python libraries.

import numpy as np
from scipy.sparse import dia_matrix # For handling diagonal matrices
from numpy.linalg import inv        # For computing the inverse of a matrix

Logistic regression

We start by describing the logistic regression model. Let be a vector of binary indicators depending on a set of covariates , for so that

independently for , with . The inverse of logit function in the above equation is equal to

In Python, we can define the inverse logit function as follow

def invlogit(eta):
    return(1.0/(1.0+np.exp(-eta)))

Iteratively Reweighted Least Square (IRLS)

Let be the design matrix and let be the vector containing the probabilities of having a positive outcome. At the -th step the IRLS algorithm iteratively updates the coefficient as follow

where and where . We proceed in this way until convergence.

In the following chunk of code, we provide a simple implementation in Python of these updating equations. We remark that this code might be highly inefficient and it serves only for didactical pourposes.

def IRLS_logit(y, X, max_iter = 500):

    # Initialization
    n,p = X.shape
    W    = dia_matrix((n,n))        # W is initialized
    beta = inv(X.T * X) * X.T * y   # OLS estimate as initialization
    
    # IRLS
    for iter in range(max_iter):
        pi = invlogit(X * beta)                # Evaluate the probabilities
        W.setdiag(np.multiply(pi,(1 - pi)).A1) # Set the diagonal
        # Updating beta
        beta_star = beta + inv(X.T * W * X) * X.T * (y - pi)
        # Check for convergence
        error = max(abs((beta_star - beta)/beta_star))
        if error < 1e-10:
            print("Convergence reached after",iter+1,"iterations")
            return({'Estimate':beta_star,'Iter':iter})
        # If the convergence criterium is not satisfied, continue
        beta = beta_star
    print("Maximum iteration reached without convergence")

Expectation Maximization

Following Scott and Sun (2013) we can obtain a different iterative procedure based on a Polya-Gamma data augmentation. At the -th step the EM algorithm iteratively updates the coefficient as follow

where . We proceed in this way until convergence. The Python code is

def EM_logit(y, X, max_iter = 500):

    # Initialization
    n,p   = X.shape
    Omega = dia_matrix((n,n))        # Omega is initialized
    beta  = inv(X.T * X) * X.T * y   # OLS estimate as initialization
    
    # EM
    for iter in range(max_iter):
        eta = (X * beta).A1                   # Evaluate the probabilities
        Omega.setdiag(np.tanh(eta/2)/(2*eta)) # Set the diagonal
        # Updating beta
        beta_star = inv(X.T * Omega * X) * X.T * (y - 0.5)
        # Check for convergence
        error = max(abs((beta_star - beta)/beta_star))
        if error < 1e-10:
            print("Convergence reached after",iter+1,"iterations")
            return({'Estimate':beta_star,'Iter':iter})
        # If the convergence criterium is not satisfied, continue
        beta = beta_star
    print("Maximum iteration reached without convergence")

Simulation study

We simulate some fictitious data to check whether our algorithms are working properly. We assume that

independently for , with . The entries of $x_i$, are sampled uniformly over . In Python, we get

n        = 50000
X        = np.matrix([np.repeat(1,n),np.random.uniform(-2,2,n)]).T
beta_sim = np.matrix([1,1]).T
pi_sim   = invlogit(X * beta_sim)
y        = np.random.binomial(1,pi_sim)

Then, we estimate the MLE via IRLS

out_IRLS = IRLS_logit(y,X)
out_IRLS['Estimate']

and via the EM algorithm

out_EM = EM_logit(y,X)
out_EM['Estimate']

In both cases we reached convergence in a reasonable time.