Likelihood, score, and fitting
Section of Logistic regression, part of Statistical learning
Likelihood, score, and fitting
| Type | Section |
| Taught in | Logistic regression |
| Part | Statistical learning |
During fitting, \beta is the candidate parameter and p_i=\sigma(\beta^\top x_i) is its model probability, not yet an estimated probability. Treat the design vectors, including their intercept entries, as fixed. Conditional independence of the responses then gives one Bernoulli factor per observation: p_i if the label is one and 1-p_i if it is zero.
Definition 1 (Binary logistic regression) The Bernoulli likelihood and log-likelihood are
\begin{align*} L(\beta)&=\prod_{i=1}^n p_i^{y_i}(1-p_i)^{1-y_i},\\ \ell(\beta)&=\sum_{i=1}^n\bigl\{y_i\log p_i+(1-y_i)\log(1-p_i)\bigr\}, \end{align*}
with p_i=\sigma(\beta^\top x_i) and y_i\in\{0,1\}.
The log-likelihood is the object to maximise. The logit link says \log(p_i/(1-p_i))=\beta^\top x_i=:\eta_i. Solving for the two log probabilities,
\log p_i=\eta_i-\log(1+e^{\eta_i}), \qquad \log(1-p_i)=-\log(1+e^{\eta_i}),
because p_i=e^{\eta_i}/(1+e^{\eta_i}) and 1-p_i=1/(1+e^{\eta_i}). Substituting those expressions removes p_i from the criterion:
\begin{align*} \ell(\beta) &=\sum_{i=1}^n\Bigl\{y_i\bigl(\eta_i-\log(1+e^{\eta_i})\bigr)+(1-y_i)\bigl(-\log(1+e^{\eta_i})\bigr)\Bigr\}\\ &=\sum_{i=1}^n\bigl\{y_i\eta_i-\log(1+e^{\eta_i})\bigr\}. \end{align*}
Score and Hessian
Main article: Score and Hessian of logistic regression
There is no closed form for the maximiser, so the next step is the derivative. Differentiate the reduced log-likelihood with \partial\eta_i/\partial\beta=x_i and
\frac{\partial}{\partial\eta_i}\log(1+e^{\eta_i})=\frac{e^{\eta_i}}{1+e^{\eta_i}}=p_i.
The score, a vector of length p+1, is
\frac{\partial\ell}{\partial\beta} =\sum_{i=1}^n(y_i-p_i)\,x_i =X^\top(y-p),
where X is the n\times(p+1) design matrix, y is the vector of labels, and p is the vector of fitted probabilities. Differentiating again, \partial p_i/\partial\beta=p_i(1-p_i)\,x_i by the logistic derivative \sigma'(\eta_i)=p_i(1-p_i). The Hessian, of size (p+1)\times(p+1), is therefore
\frac{\partial^2\ell}{\partial\beta\,\partial\beta^\top} =-\sum_{i=1}^n p_i(1-p_i)\,x_i x_i^\top =-X^\top W X,
where W=\mathrm{diag}\bigl(p_i(1-p_i)\bigr). The identity p_i(1-p_i)=\sigma'(\eta_i) is what puts the weights on the diagonal of W.
Theorem 1 (A finite stationary point is a global maximiser) At every finite coefficient vector the weights p_i(1-p_i) are positive. The Hessian -X^\top W X is negative semi-definite, so \ell is concave. If X has full column rank, W is positive definite on the column space of X and the Hessian is negative definite, so \ell is strictly concave. Any finite stationary point is then the unique global maximiser.
Concavity does not prove that a finite stationary point exists, and it does not prove that an unguarded numerical iteration will find one.
Complete separation is a concrete failure. Suppose there is a direction v such that v^\top x_i>0 whenever y_i=1 and v^\top x_i<0 whenever y_i=0. Move the coefficient vector farther along that direction: the probabilities approach the observed zeroes and ones, and the log-likelihood approaches its upper bound 0 without attaining it at a finite vector. The bound holds because each Bernoulli factor is at most 1, so L(\beta)\le 1 and \ell(\beta)\le 0, with equality only in the limit where every fitted probability is exactly the observed label. The likelihood stays bounded. The coefficients diverge. Quasi-complete separation can also prevent a finite maximum likelihood estimate. Check convergence and separation before interpreting standard errors.
Iteratively reweighted least squares
Main article: Iteratively reweighted least squares
Newton–Raphson updates
\beta\leftarrow\beta+(X^\top W X)^{-1}X^\top(y-p).
The update is the normal equation of a weighted least-squares problem, so each iteration can be computed by a regression routine. Define the working response
z=X\beta+W^{-1}(y-p).
Then
\begin{align*} X^\top W z &=X^\top W X\,\beta+X^\top(y-p),\\ (X^\top W X)^{-1}X^\top W z &=\beta+(X^\top W X)^{-1}X^\top(y-p), \end{align*}
which is exactly the Newton step. Each step is a weighted least-squares fit of z on X, with weights p_i(1-p_i). The weights are functions of the current \beta, through p_i, and they are recomputed every iteration. This is iteratively reweighted least squares. It typically converges in a handful of iterations.
For logistic regression, Fisher scoring and Newton–Raphson coincide, because the Hessian -X^\top W X does not depend on y except through the current probabilities, and those probabilities are functions of \beta and X. The expected Hessian and the observed Hessian are the same matrix once the variance of a Bernoulli trial is taken to be p_i(1-p_i). The working response is a device that makes one Newton step look like weighted least squares. It is recomputed from the current fit, and it is not itself the binary observation.
Example 1 (Newton’s method on four points) Question. Fit \log(p/(1-p))=\beta_0+\beta_1 x to x=(0,1,2,3) and y=(0,1,0,1).
Solution. The design matrix has rows (1,x_i). Start at \beta=0: every p_i=1/2, W=(1/4)I, and y-p=(-1/2,1/2,-1/2,1/2). The score is
X^\top(y-p) =\begin{pmatrix} -1/2+1/2-1/2+1/2\\ 0\cdot(-1/2)+1\cdot(1/2)+2\cdot(-1/2)+3\cdot(1/2) \end{pmatrix} =\begin{pmatrix}0\\1\end{pmatrix}.
Also X^\top X=\begin{pmatrix}4&6\\6&14\end{pmatrix}, so
X^\top W X=\frac14 X^\top X=\begin{pmatrix}1&1.5\\1.5&3.5\end{pmatrix}.
The determinant is 3.5-2.25=1.25, and
(X^\top W X)^{-1} =\frac{1}{1.25}\begin{pmatrix}3.5&-1.5\\-1.5&1\end{pmatrix}.
The Newton step is that inverse times (0,1)^\top,
\Delta\beta=\frac{1}{1.25}\begin{pmatrix}-1.5\\1\end{pmatrix} =\begin{pmatrix}-1.2\\0.8\end{pmatrix}.
After this first step \beta=(-1.2,0.8)^\top. The log-likelihood has risen from 4\log(1/2)=-2.7726 to -2.3526. Two more steps give (-1.356,0.904)^\top and (-1.3623,0.9082)^\top, with \ell=-2.3475, after which the score is zero to four decimal places. The asymptotic standard errors, the square roots of the diagonal of (X^\top W X)^{-1} at the solution, are 1.97 and 1.09. With four observations, nothing is distinguishable from zero at the usual scale.
Used in
See also
Binary logistic regression · A finite stationary point is a global maximiser · From a linear model to class probabilities · Logistic and probit models · Score and Hessian of logistic regression · Iteratively reweighted least squares · Interpretation and the decision boundary · Linear discriminant analysis implies a logistic model