Poisson Regression › Formula

The foundation of Poisson regression lies in the Poisson distribution, which describes the probability of observing exactly $k$ events when events occur independently at a constant average rate $\lambda$.

The probability mass function (PMF) of the Poisson distribution, which forms the basis of Poisson regression, is given by:

$$P(Y = k \mid \lambda) = \frac{e^{-\lambda} \lambda^k}{k!}$$

where:

• $Y$ is a discrete random variable representing the count of events (the outcome variable we observe)
• $k$ is a non-negative integer ($k = 0, 1, 2, \ldots$) representing the observed count (the specific value taken by $Y$)
• $\lambda$ is the rate parameter (the expected number of events in a given interval, where $\lambda > 0$; also equal to both the mean and variance of the distribution)
• $k!$ denotes the factorial of $k$ (i.e., $k! = k \times (k-1) \times \cdots \times 2 \times 1$, with $0! = 1$ by definition)
• $e$ is Euler's number (approximately 2.71828; the base of the natural logarithm)

This formula expresses the probability of observing exactly $k$ events in a fixed interval of time or space, assuming events occur independently and at a constant average rate $\lambda$.

In the context of Poisson regression, $\lambda$ is not a fixed value but is modeled as a function of predictor variables, allowing us to understand how different factors influence the expected event count. This distribution naturally handles count data because it only assigns positive probability to non-negative integers.

In Poisson regression, our goal is to model the relationship between a set of predictor variables and the expected rate at which events occur, denoted by the parameter $\lambda$. This parameter represents the mean number of occurrences for a given set of predictors. Because $\lambda$ needs to be positive (since you cannot have a negative rate of occurrence), we use a logarithmic link function to relate the linear combination of predictors to $\lambda$. This is known as the log-link function, and it is a defining feature of Poisson regression.

The model is specified as follows:

$$\log(\lambda_i) = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + \cdots + \beta_p x_{ip}$$

where:

• $\lambda_i$ is the expected count (rate parameter) for observation $i$ (where $\lambda_i > 0$; represents the mean of the Poisson distribution for observation $i$)
• $\beta_0$ is the intercept (the log of the expected count when all predictors equal zero)
• $\beta_1, \ldots, \beta_p$ are the coefficients for each predictor (representing the change in log-count per unit change in the predictor, holding all other predictors constant)
• $x_{i1}, \ldots, x_{ip}$ are the values of the predictors for observation $i$ (the feature values for the $i$-th observation)
• $p$ is the number of predictor variables in the model (total count of features excluding the intercept)
• $i$ indexes the observations from $1$ to $n$, where $n$ is the total number of observations

This logarithmic transformation guarantees that the predicted rate $\lambda_i$ is positive, regardless of the values of the predictors or coefficients. By exponentiating both sides, we see that:

$$\lambda_i = \exp(\beta_0 + \beta_1 x_{i1} + \cdots + \beta_p x_{ip})$$

where:

• $\exp(\cdot)$ is the exponential function (the inverse of the natural logarithm; maps any real number to a positive value)
• $\lambda_i$ is guaranteed to be positive for any values of $\beta_0, \beta_1, \ldots, \beta_p$ and $x_{i1}, \ldots, x_{ip}$

The exponential function maps any real number to a positive value, which is exactly what we need for modeling count data. This also means that the effect of each predictor is multiplicative on the expected count: a one-unit increase in $x_{ij}$ multiplies $\lambda_i$ by $e^{\beta_j}$, holding other variables constant.

For computational efficiency and clarity, we often rewrite the model in matrix notation. Let $\mathbf{X}$ be the design matrix containing all predictor variables (with a column of ones for the intercept), and let $\boldsymbol{\beta}$ be the vector of coefficients. Then, for all $n$ observations:

$$\log(\boldsymbol{\lambda}) = \mathbf{X}\boldsymbol{\beta}$$

where:

• $\boldsymbol{\lambda}$ is the $n \times 1$ vector of rate parameters for all observations: $\boldsymbol{\lambda} = [\lambda_1, \lambda_2, \ldots, \lambda_n]^T$ (the expected counts for each observation)
• $\mathbf{X}$ is the $n \times (p+1)$ design matrix where each row represents one observation and each column represents a predictor (the first column contains all ones for the intercept term)
• $\boldsymbol{\beta}$ is the $(p+1) \times 1$ vector of coefficients: $\boldsymbol{\beta} = [\beta_0, \beta_1, \ldots, \beta_p]^T$ (includes the intercept and all predictor coefficients)
• $n$ is the number of observations in the dataset (total number of data points)
• $\log(\boldsymbol{\lambda})$ denotes element-wise application of the logarithm to each element of $\boldsymbol{\lambda}$

This matrix formulation allows us to efficiently compute predictions for all data points and to fit the model using numerical optimization algorithms.

The Poisson regression model is typically estimated using maximum likelihood estimation (MLE). The likelihood function for the entire dataset is the product of the individual Poisson probabilities for each observation. For $n$ independent observations, the likelihood is:

$$L(\boldsymbol{\beta}) = \prod_{i=1}^{n} \frac{e^{-\lambda_i} \lambda_i^{y_i}}{y_i!}$$

where:

• $L(\boldsymbol{\beta})$ is the likelihood function (the probability of observing the entire dataset given the parameter vector $\boldsymbol{\beta}$)
• $y_i$ is the observed count for observation $i$ (a non-negative integer; the actual data value we observed)
• $\lambda_i$ is the expected count for observation $i$, as determined by the predictors and coefficients: $\lambda_i = \exp(\mathbf{x}_i^T \boldsymbol{\beta})$
• $\mathbf{x}_i$ is the $(p+1) \times 1$ feature vector for observation $i$ (includes 1 for the intercept and the $p$ predictor values)
• $\boldsymbol{\beta}$ is the $(p+1) \times 1$ vector of coefficients we are estimating: $\boldsymbol{\beta} = [\beta_0, \beta_1, \ldots, \beta_p]^T$
• $\prod_{i=1}^{n}$ denotes the product over all $n$ observations (assuming observations are independent)

Because products of many small probabilities can be numerically unstable, we usually work with the log-likelihood, which is the sum of the logarithms of the individual probabilities:

$$\ell(\boldsymbol{\beta}) = \sum_{i=1}^{n} \left[ -\lambda_i + y_i \log(\lambda_i) - \log(y_i!) \right]$$

where:

• $\ell(\boldsymbol{\beta})$ is the log-likelihood function (the natural logarithm of the likelihood function)
• $\sum_{i=1}^{n}$ denotes the sum over all $n$ observations (converting the product in the likelihood to a sum)
• $\log(\cdot)$ is the natural logarithm (base $e$; the inverse of the exponential function)
• $-\lambda_i$ comes from $\log(e^{-\lambda_i}) = -\lambda_i$
• $y_i \log(\lambda_i)$ comes from $\log(\lambda_i^{y_i}) = y_i \log(\lambda_i)$
• $\log(y_i!)$ comes from the factorial term in the denominator

The term $\log(y_i!)$ does not depend on the model parameters $\boldsymbol{\beta}$, so it can be omitted during the optimization process (it acts as a constant additive term that doesn't affect which $\boldsymbol{\beta}$ maximizes the function). This simplifies the log-likelihood to:

$$\ell(\boldsymbol{\beta}) = \sum_{i=1}^{n} \left[ -\lambda_i + y_i \log(\lambda_i) \right]$$

To find the values of the coefficients $\boldsymbol{\beta}$ that best fit the data, we maximize the log-likelihood function. This involves taking the partial derivatives of the log-likelihood with respect to each coefficient and setting them equal to zero.

Let's expand the steps for computing the derivative of the log-likelihood with respect to $\beta_j$.

First, recall the log-likelihood:

$$\ell(\boldsymbol{\beta}) = \sum_{i=1}^n \left[ -\lambda_i + y_i \log(\lambda_i) \right]$$

where:

• $\lambda_i = \exp(\mathbf{x}_i^T \boldsymbol{\beta})$ is the expected count for observation $i$ (a function of the coefficients we're optimizing)
• $\mathbf{x}_i^T \boldsymbol{\beta} = \sum_{j=0}^p x_{ij} \beta_j$ is the linear predictor for observation $i$ (the sum of all coefficients times their corresponding predictor values; the argument to the exponential function)
• $\mathbf{x}_i = [1, x_{i1}, x_{i2}, \ldots, x_{ip}]^T$ is the $(p+1) \times 1$ feature vector for observation $i$ (the first element is 1 for the intercept term)
• $j$ indexes the coefficients from $0$ to $p$, where $j=0$ corresponds to the intercept $\beta_0$

Take the partial derivative with respect to $\beta_j$:

$$\begin{aligned} \frac{\partial \ell}{\partial \beta_j} &= \sum_{i=1}^n \left[ -\frac{\partial \lambda_i}{\partial \beta_j} + y_i \frac{1}{\lambda_i} \frac{\partial \lambda_i}{\partial \beta_j} \right] \end{aligned}$$

where:

• $\frac{\partial \ell}{\partial \beta_j}$ is the partial derivative of the log-likelihood with respect to coefficient $\beta_j$ (the gradient component for the $j$-th parameter)
• The first term $-\frac{\partial \lambda_i}{\partial \beta_j}$ comes from differentiating $-\lambda_i$
• The second term $y_i \frac{1}{\lambda_i} \frac{\partial \lambda_i}{\partial \beta_j}$ comes from differentiating $y_i \log(\lambda_i)$ using the chain rule

Now, compute $\frac{\partial \lambda_i}{\partial \beta_j}$ using the chain rule:

$$\begin{aligned} \lambda_i &= \exp(\mathbf{x}_i^T \boldsymbol{\beta}) \\ \frac{\partial \lambda_i}{\partial \beta_j} &= \exp(\mathbf{x}_i^T \boldsymbol{\beta}) \cdot \frac{\partial (\mathbf{x}_i^T \boldsymbol{\beta})}{\partial \beta_j} \\ &= \exp(\mathbf{x}_i^T \boldsymbol{\beta}) \cdot x_{ij} \\ &= \lambda_i x_{ij} \end{aligned}$$

where:

• $\frac{\partial \lambda_i}{\partial \beta_j}$ is the partial derivative of the expected count with respect to coefficient $\beta_j$ (how $\lambda_i$ changes when $\beta_j$ changes)
• The first step uses the chain rule: $\frac{d}{dx}e^{f(x)} = e^{f(x)} \cdot f'(x)$
• $\frac{\partial (\mathbf{x}_i^T \boldsymbol{\beta})}{\partial \beta_j} = x_{ij}$ because $\mathbf{x}_i^T \boldsymbol{\beta} = \sum_{k=0}^p x_{ik} \beta_k$, and only the $j$-th term depends on $\beta_j$
• $x_{ij}$ is the $j$-th element of the feature vector $\mathbf{x}_i$ (the value of predictor $j$ for observation $i$)
• For $j=0$, $x_{i0} = 1$ (the intercept term)
• For $j > 0$, $x_{ij}$ is the value of the $j$-th predictor variable for observation $i$

Substitute this result back into the derivative:

$$\begin{aligned} \frac{\partial \ell}{\partial \beta_j} &= \sum_{i=1}^n \left[ -\lambda_i x_{ij} + y_i \frac{1}{\lambda_i} \lambda_i x_{ij} \right] \\ &= \sum_{i=1}^n \left[ -\lambda_i x_{ij} + y_i x_{ij} \right] \\ &= \sum_{i=1}^n x_{ij} (y_i - \lambda_i) \end{aligned}$$

where:

• In the first line, we substitute $\frac{\partial \lambda_i}{\partial \beta_j} = \lambda_i x_{ij}$ into both terms
• In the second line, we simplify $y_i \frac{1}{\lambda_i} \lambda_i x_{ij} = y_i x_{ij}$ (the $\lambda_i$ terms cancel)
• In the third line, we factor out $x_{ij}$ from both terms to get the final form

To find the maximum likelihood estimates, set the derivative equal to zero and solve for $\beta_j$:

$$\sum_{i=1}^n x_{ij} (y_i - \lambda_i) = 0 \quad \text{for all } j = 0, 1, \ldots, p$$

where:

• This equation must hold for each coefficient $\beta_j$ (including the intercept $\beta_0$ when $j=0$)
• $y_i - \lambda_i$ represents the residual (observed count minus expected count) for observation $i$ (the prediction error)
• The equation states that the weighted sum of residuals must equal zero, where weights are the predictor values $x_{ij}$
• This is the first-order condition for a maximum (the gradient equals zero)

This yields a system of $(p+1)$ equations (one for each coefficient, including the intercept), which are solved numerically to obtain the maximum likelihood estimates $\hat{\beta}_0, \hat{\beta}_1, \ldots, \hat{\beta}_p$.

These equations are nonlinear in $\boldsymbol{\beta}$ and generally do not have a closed-form solution. Therefore, we solve them numerically using iterative optimization algorithms. The most common approaches are:

1. Newton-Raphson method: Uses both first and second derivatives (the Hessian matrix) to find the optimal coefficients
2. Iteratively reweighted least squares (IRLS): Recasts each iteration as a weighted least squares problem

The numerical solution process involves:

• Starting with initial guesses for the coefficients (e.g., $\boldsymbol{\beta}^{(0)} = \mathbf{0}$, where $\mathbf{0}$ is a vector of zeros)
• Computing predicted values $\lambda_i$ using the current coefficient estimates: $\lambda_i^{(t)} = \exp(\mathbf{x}_i^T \boldsymbol{\beta}^{(t)})$, where $t$ denotes the iteration number
• Updating coefficients using gradient information (the partial derivatives): $\boldsymbol{\beta}^{(t+1)} = \boldsymbol{\beta}^{(t)} + \Delta\boldsymbol{\beta}$, where $\Delta\boldsymbol{\beta}$ is the update step computed from the gradient
• Repeating until convergence (when $\|\boldsymbol{\beta}^{(t+1)} - \boldsymbol{\beta}^{(t)}\| < \epsilon$ for some small tolerance $\epsilon$, typically $\epsilon = 10^{-4}$ or $10^{-6}$)

When the algorithm converges, we obtain the maximum likelihood estimates $\hat{\beta}_0$, $\hat{\beta}_1$, $\ldots$, $\hat{\beta}_p$. These are the values that make the observed data most probable under the Poisson regression model.

The key insight is that the system of equations $\sum_{i=1}^{n} x_{ij} (y_i - \lambda_i) = 0$ represents the condition that the residuals $(y_i - \lambda_i)$ are orthogonal to each predictor variable, weighted by the predictor values. This ensures that the model captures the relationship between predictors and outcomes in a way that maximizes the likelihood of observing the actual data.

Note: The use of the log-link function and the exponential mean structure distinguishes Poisson regression from ordinary linear regression. This approach ensures that predictions are non-negative and that the model is well-suited for count data, where the variance typically increases with the mean.

To better understand how maximum likelihood estimation works in practice, let's visualize the iterative optimization process. The following plots show how the algorithm searches for the optimal coefficients by iteratively improving the log-likelihood.

The optimization process demonstrates several key insights about maximum likelihood estimation in Poisson regression. First, the log-likelihood surface is smooth and unimodal (has a single peak), which ensures that iterative algorithms reliably converge to the global maximum regardless of the starting point. This is a desirable property that makes Poisson regression computationally stable.

Second, the algorithm makes rapid progress in early iterations when far from the optimum, then slows down as it approaches the maximum. This is characteristic of gradient-based methods like Newton-Raphson and IRLS, which use information about both the direction and curvature of the log-likelihood surface to move efficiently toward the solution.

Finally, convergence is achieved when the gradient (rate of change) of the log-likelihood becomes negligibly small, indicating that we've reached a point where no small adjustment to the coefficients will improve the fit. In practice, algorithms stop when the change in coefficients or log-likelihood between iterations falls below a specified tolerance threshold.

The Poisson regression model has several important mathematical properties that make it well-suited for count data.

Equidispersion property: For a Poisson distribution, the mean and variance are both equal to $\lambda$:

$$E[Y] = \text{Var}(Y) = \lambda$$

where:

• $E[Y]$ is the expected value (mean) of the count variable $Y$
• $\text{Var}(Y)$ is the variance of the count variable $Y$
• $\lambda$ is the rate parameter (constrained to be positive)

This means that as the expected count increases, so does the variability.a realistic assumption for many count processes.

Scale invariance: The model is scale-invariant under the log-link: multiplying a predictor by a constant $c$ simply shifts the corresponding coefficient by $\log(c)$. If $x_j^* = c \cdot x_j$, then $\beta_j^* = \beta_j / c$ such that $\beta_j^* x_j^* = \beta_j x_j$.

Coefficient interpretation: The coefficients have a natural interpretation as log-rate ratios. When we exponentiate a coefficient $\beta_j$, we get the rate ratio:

$$\text{Rate Ratio} = e^{\beta_j}$$

This represents the multiplicative effect on the expected count for a one-unit increase in $x_j$, holding all other variables constant. This multiplicative structure is often more appropriate for count data than the additive structure of linear regression, as it naturally handles the non-negative constraint and allows for proportional changes in rates.