Poisson Regression¶
1. Motivation¶
Poisson regression is used when the response variable is a count:
- number of clicks,
- number of accidents,
- number of emails,
- number of customers,
- number of calls,
- number of failures,
- number of events in a fixed time interval.
The response variable satisfies
$$
Y_i \in \{0,1,2,\dots\}.
$$
Linear regression is not appropriate because it may predict negative values.
Poisson regression solves this by modeling the mean as a positive quantity.
2. Poisson Distribution¶
A random variable $Y$ has Poisson distribution with parameter $\lambda>0$ if
$$
P(Y=y)=\frac{e^{-\lambda}\lambda^y}{y!},
\qquad y=0,1,2,\dots
$$
We write
$$
Y\sim \operatorname{Poisson}(\lambda).
$$
The mean and variance are
$$
\mathbb{E}[Y]=\lambda,
$$
and
$$
\operatorname{Var}(Y)=\lambda.
$$
Thus, for a Poisson random variable,
$$
\mathbb{E}[Y]=\operatorname{Var}(Y).
$$
This is an important property.
3. The Poisson Regression Model¶
Suppose we observe independent data
$$
(y_i,x_i), \qquad i=1,\dots,n,
$$
where
$$
x_i=
\begin{pmatrix}
1\\
x_{i1}\\
\vdots\\
x_{ip}
\end{pmatrix}
\in \mathbb{R}^{p+1}.
$$
We assume
$$
Y_i\sim \operatorname{Poisson}(\lambda_i),
$$
where the mean depends on the predictors.
The problem is that we need
$$
\lambda_i>0.
$$
Therefore, instead of modeling $\lambda_i$ directly as a linear function, we model its logarithm:
$$
\log(\lambda_i)=x_i^T\beta.
$$
Equivalently,
$$
\lambda_i=e^{x_i^T\beta}.
$$
This automatically guarantees
$$
\lambda_i>0.
$$
4. Why Use the Log Link?¶
If we used
$$
\lambda_i=x_i^T\beta,
$$
then the model could produce negative predicted counts.
That is impossible because
$$
\lambda_i=\mathbb{E}[Y_i]>0.
$$
The log link solves this:
$$
\log(\lambda_i)=x_i^T\beta.
$$
Since the exponential function is always positive,
$$
\lambda_i=e^{x_i^T\beta}>0.
$$
Thus Poisson regression is a generalized linear model with:
| Component |
Poisson Regression |
| Response distribution |
Poisson |
| Mean |
$\lambda_i$ |
| Link function |
$\log(\lambda_i)$ |
| Linear predictor |
$x_i^T\beta$ |
5. Interpretation of the Coefficients¶
Consider the model
$$
\log(\lambda_i)=\beta_0+\beta_1x_i.
$$
Then
$$
\lambda_i=e^{\beta_0+\beta_1x_i}.
$$
Now increase $x_i$ by one unit. The new mean is
$$
\lambda(x_i+1)=e^{\beta_0+\beta_1(x_i+1)}.
$$
Therefore,
$$
\lambda(x_i+1)
=
e^{\beta_0+\beta_1x_i+\beta_1}
=
e^{\beta_1}e^{\beta_0+\beta_1x_i}.
$$
So
$$
\lambda(x_i+1)=e^{\beta_1}\lambda(x_i).
$$
Hence,
$$
\frac{\lambda(x_i+1)}{\lambda(x_i)}=e^{\beta_1}.
$$
Therefore, a one-unit increase in $x_i$ multiplies the expected count by
$$
e^{\beta_1}.
$$
So:
- if $\beta_1>0$, the expected count increases;
- if $\beta_1<0$, the expected count decreases;
- if $\beta_1=0$, the predictor has no effect.
For example, if
$$
\beta_1=0.2,
$$
then
$$
e^{0.2}\approx 1.221.
$$
So a one-unit increase in $x$ increases the expected count by approximately $22.1\%$.
6. Likelihood Function¶
Since
$$
Y_i\sim \operatorname{Poisson}(\lambda_i),
$$
we have
$$
P(Y_i=y_i)
=
\frac{e^{-\lambda_i}\lambda_i^{y_i}}{y_i!}.
$$
Assuming independence, the likelihood function is
$$
L(\beta)
=
\prod_{i=1}^n
\frac{e^{-\lambda_i}\lambda_i^{y_i}}{y_i!}.
$$
Since
$$
\lambda_i=e^{x_i^T\beta},
$$
we get
$$
L(\beta)
=
\prod_{i=1}^n
\frac{
e^{-e^{x_i^T\beta}}
\left(e^{x_i^T\beta}\right)^{y_i}
}{y_i!}.
$$
Using
$$
\left(e^{x_i^T\beta}\right)^{y_i}
=
e^{y_i x_i^T\beta},
$$
we obtain
$$
L(\beta)
=
\prod_{i=1}^n
\frac{
e^{-e^{x_i^T\beta}}
e^{y_i x_i^T\beta}
}{y_i!}.
$$
7. Log-Likelihood Function¶
Taking logarithms,
$$
\ell(\beta)=\log L(\beta).
$$
Therefore,
$$
\ell(\beta)
=
\sum_{i=1}^n
\log
\left(
\frac{e^{-\lambda_i}\lambda_i^{y_i}}{y_i!}
\right).
$$
So
$$
\ell(\beta)
=
\sum_{i=1}^n
\left(
-\lambda_i+y_i\log(\lambda_i)-\log(y_i!)
\right).
$$
Since
$$
\lambda_i=e^{x_i^T\beta}
$$
and
$$
\log(\lambda_i)=x_i^T\beta,
$$
we get
$$
\ell(\beta)
=
\sum_{i=1}^n
\left(
-e^{x_i^T\beta}
+
y_i x_i^T\beta
-
\log(y_i!)
\right).
$$
Thus,
$$
\boxed{
\ell(\beta)
=
\sum_{i=1}^n
\left[
y_i x_i^T\beta
-
e^{x_i^T\beta}
-
\log(y_i!)
\right].
}
$$
Since $\log(y_i!)$ does not depend on $\beta$, for optimization we may ignore it.
Thus we maximize
$$
\tilde{\ell}(\beta)
=
\sum_{i=1}^n
\left[
y_i x_i^T\beta
-
e^{x_i^T\beta}
\right].
$$
8. Score Function¶
The score function is the gradient of the log-likelihood:
$$
U(\beta)=\nabla_\beta \ell(\beta).
$$
We have
$$
\ell(\beta)
=
\sum_{i=1}^n
\left[
y_i x_i^T\beta
-
e^{x_i^T\beta}
-
\log(y_i!)
\right].
$$
Differentiate term by term.
First,
$$
\nabla_\beta (y_i x_i^T\beta)=y_i x_i.
$$
Second,
$$
\nabla_\beta e^{x_i^T\beta}
=
e^{x_i^T\beta}x_i.
$$
Therefore,
$$
\nabla_\beta
\left[
-y_i?
\right]
$$
More precisely,
$$
\nabla_\beta
\left[
y_i x_i^T\beta
-
e^{x_i^T\beta}
\right]
=
y_i x_i
-
e^{x_i^T\beta}x_i.
$$
Hence,
$$
U(\beta)
=
\sum_{i=1}^n
x_i
\left(
y_i-e^{x_i^T\beta}
\right).
$$
Since
$$
\lambda_i=e^{x_i^T\beta},
$$
we get
$$
\boxed{
U(\beta)
=
\sum_{i=1}^n
x_i(y_i-\lambda_i).
}
$$
In matrix form, let
$$
y=
\begin{pmatrix}
y_1\\
\vdots\\
y_n
\end{pmatrix},
\qquad
\lambda=
\begin{pmatrix}
\lambda_1\\
\vdots\\
\lambda_n
\end{pmatrix}.
$$
Then
$$
\boxed{
U(\beta)=X^T(y-\lambda).
}
$$
The maximum likelihood estimator satisfies
$$
X^T(y-\lambda)=0.
$$
That is,
$$
X^T y = X^T \lambda.
$$
But since
$$
\lambda_i=e^{x_i^T\beta},
$$
this is a nonlinear system.
Therefore, there is generally no closed-form solution.
9. Hessian Matrix¶
The Hessian matrix is
$$
H(\beta)=\nabla_\beta^2 \ell(\beta).
$$
From the score function,
$$
U(\beta)=\sum_{i=1}^n x_i(y_i-\lambda_i).
$$
Since $y_i$ is constant with respect to $\beta$,
$$
\nabla_\beta(y_i-\lambda_i)
=
-\nabla_\beta \lambda_i.
$$
But
$$
\lambda_i=e^{x_i^T\beta},
$$
so
$$
\nabla_\beta \lambda_i
=
\lambda_i x_i.
$$
Therefore,
$$
\nabla_\beta
\left[
x_i(y_i-\lambda_i)
\right]
=
-\lambda_i x_i x_i^T.
$$
Hence,
$$
\boxed{
H(\beta)
=
-\sum_{i=1}^n
\lambda_i x_i x_i^T.
}
$$
In matrix form,
$$
\boxed{
H(\beta)
=
-X^T W X,
}
$$
where
$$
W=
\operatorname{diag}(\lambda_1,\dots,\lambda_n).
$$
Because $W$ is positive diagonal, $X^TWX$ is positive semidefinite.
Therefore,
$$
H(\beta)
$$
is negative semidefinite.
So the log-likelihood is concave.
This is good: any local maximum is a global maximum.
The Fisher information matrix is
$$
I(\beta)
=
-\mathbb{E}[H(\beta)].
$$
For Poisson regression,
$$
H(\beta)=-X^T W X,
$$
where
$$
W=\operatorname{diag}(\lambda_1,\dots,\lambda_n).
$$
Since $\lambda_i$ is determined by $x_i$ and $\beta$, we get
$$
\boxed{
I(\beta)=X^T W X.
}
$$
This matrix is used to approximate the variance of the MLE:
$$
\widehat{\operatorname{Var}}(\hat{\beta})
=
(X^T\widehat{W}X)^{-1},
$$
where
$$
\widehat{W}
=
\operatorname{diag}(\hat{\lambda}_1,\dots,\hat{\lambda}_n).
$$
Therefore,
$$
\operatorname{SE}(\hat{\beta}_j)
=
\sqrt{
\left[
(X^T\widehat{W}X)^{-1}
\right]_{jj}
}.
$$
11. Newton-Raphson Method¶
To maximize $\ell(\beta)$, Newton's method uses
$$
\beta^{(new)}
=
\beta^{(old)}
-
H(\beta^{(old)})^{-1}U(\beta^{(old)}).
$$
For Poisson regression,
$$
U(\beta)=X^T(y-\lambda),
$$
and
$$
H(\beta)=-X^TWX.
$$
Therefore,
$$
\beta^{(new)}
=
\beta^{(old)}
+
(X^TWX)^{-1}X^T(y-\lambda).
$$
Thus,
$$
\boxed{
\beta^{(new)}
=
\beta^{(old)}
+
(X^TWX)^{-1}X^T(y-\lambda).
}
$$
This is the Newton update for Poisson regression.
