Generalized Linear Models › Example: Logistic Regression

Let's work through a detailed example of logistic regression, which is one of the most common GLMs. Suppose we want to predict whether a student will pass an exam based on their study hours and previous test scores.

Data:

• Study hours: [2, 3, 4, 5, 6, 7, 8]
• Previous scores: [60, 65, 70, 75, 80, 85, 90]
• Pass (1) or Fail (0): [0, 0, 1, 1, 1, 1, 1]

This is our model specification:

• Distribution: Binomial
• Link function: Logit
• Model: $\log\left(\frac{p}{1-p}\right) = \beta_0 + \beta_1 \cdot \text{hours} + \beta_2 \cdot \text{scores}$

Let's solve this step by step.

Step 1: Set up the design matrix

First, we organize our data into a matrix form suitable for regression. Each row of the matrix corresponds to one observation (student), and each column is a feature (including a column of 1s for the intercept). For our example, the features are:

• Intercept (always 1)
• Study hours
• Previous test score

So, for our 7 students, the design matrix $\mathbf{X}$ looks like:

$$\mathbf{X} = \begin{bmatrix} 1 & 2 & 60 \\\\ 1 & 3 & 65 \\\\ 1 & 4 & 70 \\\\ 1 & 5 & 75 \\\\ 1 & 6 & 80 \\\\ 1 & 7 & 85 \\\\ 1 & 8 & 90 \end{bmatrix}$$

The corresponding response vector (whether the student passed) is:

$$\mathbf{y} = \begin{bmatrix} 0 \\\\ 0 \\\\ 1 \\\\ 1 \\\\ 1 \\\\ 1 \\\\ 1 \end{bmatrix}$$

Step 2: Initialize coefficients

We start the iterative process by choosing initial values for the coefficients (the $\beta$ parameters). A common choice is to set all coefficients to zero:

$$\boldsymbol{\beta}^{(0)} = \begin{bmatrix} 0 \\\\ 0 \\\\ 0 \end{bmatrix}$$

This means our initial model predicts the same probability for all students, regardless of their features.

Step 3: IRLS (Iteratively Reweighted Least Squares) Iteration

The IRLS algorithm updates the coefficients in a series of steps, using the current estimate to improve the next one. For each iteration $k$:

Note: The following steps are a detailed explanation of the IRLS algorithm to demonstrate the process. We will not go into the very specific details of the algorithm in this section. In reality, you would not compute the IRLS algorithm by hand. Modern statistical and machine learning libraries, including NumPy and scikit-learn as well as R packages, implement this algorithm efficiently.

1. Calculate the linear predictor:
   Compute $\boldsymbol{\eta}^{(k)} = \mathbf{X}\boldsymbol{\beta}^{(k)}$.

   where:

   • $\boldsymbol{\eta}^{(k)}$ is the vector of linear predictors at iteration $k$ (one entry per observation)
   • $\mathbf{X}$ is the $n \times (p+1)$ design matrix (each row is one observation)
   • $\boldsymbol{\beta}^{(k)}$ is the coefficient vector at iteration $k$
   • For the first iteration with $\boldsymbol{\beta}^{(0)} = \mathbf{0}$, all entries in $\boldsymbol{\eta}^{(0)}$ equal 0

2. Calculate fitted values (predicted probabilities):
   Apply the inverse link function (the logistic sigmoid for logistic regression) to get the predicted probability of passing for each student:

   $$\boldsymbol{\mu}^{(k)} = g^{-1}(\boldsymbol{\eta}^{(k)}) = \frac{1}{1 + e^{-\boldsymbol{\eta}^{(k)}}}$$

   where:

   • $\boldsymbol{\mu}^{(k)}$ is the vector of fitted values (predicted probabilities) at iteration $k$
   • $g^{-1}(\cdot)$ is the inverse link function (for logistic regression, this is the logistic sigmoid function)
   • The formula is applied element-wise: $\mu_i^{(k)} = \frac{1}{1 + e^{-\eta_i^{(k)}}}$ for each observation $i$
   • For the first iteration with $\boldsymbol{\eta}^{(0)} = \mathbf{0}$, all predicted probabilities equal 0.5 (since $\frac{1}{1+e^0} = \frac{1}{2}$)

3. Calculate the working response:
   This step creates a "pseudo-response" that allows us to use weighted least squares:

   $$\mathbf{z}^{(k)} = \boldsymbol{\eta}^{(k)} + \frac{\mathbf{y} - \boldsymbol{\mu}^{(k)}}{\boldsymbol{\mu}^{(k)}(1 - \boldsymbol{\mu}^{(k)})}$$

   where:

   • $\mathbf{z}^{(k)}$ is the vector of working responses at iteration $k$ (adjusted target values for the weighted least squares step)
   • $\mathbf{y}$ is the vector of observed responses (actual outcomes: 0 or 1 for each observation)
   • $\mathbf{y} - \boldsymbol{\mu}^{(k)}$ is the residual vector (difference between actual and predicted values)
   • $\boldsymbol{\mu}^{(k)}(1 - \boldsymbol{\mu}^{(k)})$ is computed element-wise and represents the variance of the Bernoulli distribution for each observation (i.e., $\text{Var}(Y_i) = \mu_i(1-\mu_i)$)
   • Division is performed element-wise

4. Calculate the weights:
   The weights for each observation are based on the variance of the predicted probability:

   $$\mathbf{W}^{(k)} = \text{diag}(\boldsymbol{\mu}^{(k)}(1 - \boldsymbol{\mu}^{(k)}))$$

   where:

   • $\mathbf{W}^{(k)}$ is an $n \times n$ diagonal weight matrix at iteration $k$
   • $\text{diag}(\mathbf{v})$ creates a diagonal matrix with the elements of vector $\mathbf{v}$ on the diagonal
   • The diagonal entries are $W_{ii}^{(k)} = \mu_i^{(k)}(1 - \mu_i^{(k)})$, which is the variance of the Bernoulli distribution for observation $i$
   • Observations with predicted probabilities near 0.5 receive higher weights (maximum variance of 0.25), while those near 0 or 1 receive lower weights

5. Update the coefficients:
   Solve the weighted least squares problem to get the new coefficient estimates:

   $$\boldsymbol{\beta}^{(k+1)} = (\mathbf{X}^T\mathbf{W}^{(k)}\mathbf{X})^{-1}\mathbf{X}^T\mathbf{W}^{(k)}\mathbf{z}^{(k)}$$

   where:

   • $\boldsymbol{\beta}^{(k+1)}$ is the updated coefficient vector for iteration $k+1$
   • $\mathbf{X}^T$ is the transpose of the design matrix
   • $(\mathbf{X}^T\mathbf{W}^{(k)}\mathbf{X})^{-1}$ is the inverse of the weighted cross-product matrix (analogous to $(\mathbf{X}^T\mathbf{X})^{-1}$ in ordinary least squares)
   • This formula is the solution to the weighted least squares problem: minimize $(\mathbf{z}^{(k)} - \mathbf{X}\boldsymbol{\beta})^T\mathbf{W}^{(k)}(\mathbf{z}^{(k)} - \mathbf{X}\boldsymbol{\beta})$

Step 4: Check for convergence

After each iteration, check if the coefficients have stabilized:

$$||\boldsymbol{\beta}^{(k+1)} - \boldsymbol{\beta}^{(k)}|| < \epsilon$$

where:

• $||\cdot||$ is the Euclidean norm (also called the L2 norm), which measures the magnitude of a vector: $||\mathbf{v}|| = \sqrt{v_1^2 + v_2^2 + \cdots + v_p^2}$
• $\boldsymbol{\beta}^{(k+1)}$ is the coefficient vector at iteration $k+1$ (the newly updated coefficients)
• $\boldsymbol{\beta}^{(k)}$ is the coefficient vector at iteration $k$ (the previous coefficients)
• $\boldsymbol{\beta}^{(k+1)} - \boldsymbol{\beta}^{(k)}$ is the change in coefficients between iterations
• $\epsilon$ is a small tolerance threshold (e.g., $\epsilon = 10^{-6}$), which defines how close to convergence we need to be

If the change is smaller than $\epsilon$, the algorithm has converged. Otherwise, repeat the steps above with the new coefficients. The process continues until the change in coefficients is very small, indicating that we've found the maximum likelihood estimates.

Summary of what happens

• We start with all coefficients at zero (no effect from features).
• At each step, we use the current model to predict probabilities, then adjust the coefficients to better fit the observed outcomes, giving more weight to observations where the model is less certain.
• This process repeats, gradually improving the fit, until the coefficients stop changing significantly.