Bernoulli Logit Score Equations

From Binary Outcomes to a Logit Likelihood

In a binary-response regression, each observation records an outcome of zero or one, while its success probability depends on covariates. The logit model turns a linear predictor into that probability. Starting with the individual Bernoulli likelihood makes the coefficient score and its curvature especially transparent: each residual is the observed outcome minus its fitted probability.

Assume \(n\) independent observations, with fixed covariate vectors \(x_i\in\mathbb{R}^p\) and responses \(y_i\in\{0,1\}\), for \(i=1,\ldots,n\). Let \(\beta\in\mathbb{R}^p\) be the coefficient vector. The linear predictor is \(\eta_i=x_i^\mathsf{T}\beta\), and the logit link sets the success probability to \(\pi_i\), where \(\operatorname{logit}(\pi_i)=\log\{\pi_i/(1-\pi_i)\}=\eta_i\). Thus \(\pi_i=\exp(\eta_i)/(1+\exp(\eta_i))\). For finite \(\beta\), each \(\pi_i\) lies strictly between zero and one.

Definition: Under the Bernoulli logit regression model, the responses are independent with \(Y_i\sim\operatorname{Bernoulli}(\pi_i)\), where \(\pi_i=\exp(x_i^\mathsf{T}\beta)/\{1+\exp(x_i^\mathsf{T}\beta)\}\). For observed responses \(y_1,\ldots,y_n\), the coefficient score is the gradient \(U(\beta)=\partial\ell(\beta)/\partial\beta\), and the Hessian is \(H(\beta)=\partial^2\ell(\beta)/\partial\beta\,\partial\beta^\mathsf{T}\).

The Bernoulli mass function gives observation \(i\) likelihood contribution \(\pi_i^{y_i}(1-\pi_i)^{1-y_i}\). Independence means the likelihood for the observed data is the product of these contributions. This is an individual-level likelihood: each \(y_i\) is a zero or a one.

Deriving the Score and Hessian

The log-likelihood can first be written in terms of \(\pi_i\). Substituting the logistic expression for \(\pi_i\) then gives a particularly simple form in terms of \(\eta_i\).

\begin{aligned} L(\beta) &=\prod_{i=1}^n \pi_i^{y_i}(1-\pi_i)^{1-y_i},\\ \ell(\beta) &=\sum_{i=1}^n\left[y_i\log(\pi_i)+(1-y_i)\log(1-\pi_i)\right]\\ &=\sum_{i=1}^n\left[y_i\eta_i-\log(1+\exp(\eta_i))\right]. \end{aligned}

For the last equality, \(\log(\pi_i)=\eta_i-\log(1+\exp(\eta_i))\) and \(\log(1-\pi_i)=-\log(1+\exp(\eta_i))\). Differentiating one summand with respect to \(\eta_i\) gives \(y_i-\pi_i\), since the derivative of \(\log(1+\exp(\eta_i))\) is \(\pi_i\). Also, \(\partial\eta_i/\partial\beta=x_i\). Applying the chain rule observation by observation yields the coefficient score.

Theorem: For the Bernoulli logit model with fixed covariates and independent responses, the log-likelihood, coefficient score, and Hessian are

\begin{aligned} \ell(\beta) &=\sum_{i=1}^n\left[y_i x_i^\mathsf{T}\beta-\log\{1+\exp(x_i^\mathsf{T}\beta)\}\right],\\ U(\beta) &=\sum_{i=1}^n x_i(y_i-\pi_i),\\ H(\beta) &=-\sum_{i=1}^n \pi_i(1-\pi_i)x_ix_i^\mathsf{T}. \end{aligned}

Here \(\pi_i=\exp(x_i^\mathsf{T}\beta)/\{1+\exp(x_i^\mathsf{T}\beta)\}\). In particular, any finite unconstrained maximum-likelihood estimate must solve \(U(\beta)=0\).

To verify the Hessian, differentiate the score. The derivative of \(\pi_i\) with respect to its linear predictor is \(\pi_i(1-\pi_i)\), and \(\partial\eta_i/\partial\beta=x_i\). Therefore, the derivative of \(x_i(y_i-\pi_i)\) with respect to \(\beta^\mathsf{T}\) is \(-\pi_i(1-\pi_i)x_ix_i^\mathsf{T}\). Summing these derivatives gives the displayed Hessian. The score equation follows from the vanishing gradient at a finite maximum in the open parameter space \(\mathbb{R}^p\). This condition does not by itself assert that such a maximum exists.

Worked Example: An Intercept-Only Model

Suppose three independent outcomes are \(y_1=1\), \(y_2=0\), and \(y_3=1\), with no covariates other than an intercept. Thus \(x_i=1\), \(\eta_i=\beta\), and all observations have fitted probability \(\pi=\exp(\beta)/(1+\exp(\beta))\).

The log-likelihood and score are

\ell(\beta)=2\beta-3\log(1+\exp(\beta)), \qquad U(\beta)=2-3\pi.

Setting the score to zero gives \(\pi=2/3\), so \(\exp(\beta)=\pi/(1-\pi)=2\) and \(\hat{\beta}=\log 2\). At this value, the fitted residuals sum to \(1-2/3+0-2/3+1-2/3=0\). The Hessian is \(-3\pi(1-\pi)\), which at the solution equals \(-3(2/3)(1/3)=-2/3\). Thus the log-likelihood curves downward there, and the score equation reproduces the observed success proportion as the fitted probability.

Score Equations as Residual Balance

Let \(X\) be the \(n\)-by-\(p\) design matrix whose \(i\)th row is \(x_i^\mathsf{T}\), and let \(y\) and \(\pi\) be the vectors of responses and fitted probabilities. The score equation has the compact form \(X^\mathsf{T}(y-\pi)=0\). Its \(j\)th coordinate says that the residuals \(y_i-\pi_i\), multiplied by covariate \(j\), sum to zero. This is a balance condition over the observed data, not a claim that every residual is zero.

Interpretation: At a finite solution of the score equations, each design-matrix column has zero inner product with the residual vector \(y-\pi\). If the model includes an intercept, its column consists of ones, so the residuals sum to zero. Other coefficient equations balance residuals against their corresponding covariates.

The condition follows directly by collecting terms in the score:

U(\beta)=\sum_{i=1}^n x_i(y_i-\pi_i) =X^\mathsf{T}(y-\pi).

This interpretation is useful because it connects a fitted model to the covariates it uses. For example, with an intercept and one numerical predictor, a solution balances the overall residual sum and the predictor-weighted residual sum. Without an intercept, there is no general requirement that the unweighted residuals sum to zero.

Worked Example: Two Covariate Balance Equations

Consider six individual outcomes. The first three have covariate row \(x_i^\mathsf{T}=(1,0)\) and responses \(1,0,0\); the next three have row \((1,1)\) and responses \(1,1,0\). The first covariate is an intercept and the second indicates the second category. At \(\beta_0=-\log 2\), \(\beta_1=\log 4\), the fitted probabilities are \(1/3\) for the first three observations and \(2/3\) for the last three.

For the intercept score, the observed outcomes sum to \(3\), while the fitted probabilities sum to \(3(1/3)+3(2/3)=3\). Hence \(U_0=3-3=0\). For the category score, only the last three rows contribute: their observed outcomes sum to \(2\), and their fitted probabilities sum to \(3(2/3)=2\). Hence \(U_1=2-2=0\). The logit for the first category is \(-\log 2\); the logit for the second is \(-\log 2+\log 4=\log 2\), matching probabilities \(1/3\) and \(2/3\).

Here each observation still contributes its own Bernoulli likelihood term. The balance calculations simply collect those individual residual contributions by covariate column.

Curvature and Concavity

The Hessian describes the local curvature of the log-likelihood. For any vector \(a\in\mathbb{R}^p\), its quadratic form is

a^\mathsf{T}H(\beta)a =-\sum_{i=1}^n \pi_i(1-\pi_i)(x_i^\mathsf{T}a)^2\leq 0.

Every weight \(\pi_i(1-\pi_i)\) is strictly positive for finite \(\beta\). Thus the Hessian is negative semidefinite, and the log-likelihood is concave. If \(X\) has full column rank, then for every nonzero \(a\), at least one \(x_i^\mathsf{T}a\) is nonzero; the quadratic form is then strictly negative. In that case the Hessian is negative definite at every finite \(\beta\), and the log-likelihood is strictly concave.

Result: The Bernoulli logit log-likelihood is concave in \(\beta\). If the design matrix has full column rank, it is strictly concave, so there can be at most one finite solution to the score equations. If such a solution exists, concavity makes it a global maximum; strict concavity makes that maximum unique.

For completeness, a twice-differentiable function on \(\mathbb{R}^p\) with negative semidefinite Hessian everywhere is concave: along any line \(\beta+ta\), its second derivative is \(a^\mathsf{T}H(\beta+ta)a\leq 0\), so its derivative along the line is nonincreasing. With a negative definite Hessian, that second derivative is strictly negative for every nonzero direction, giving strict concavity along every nonconstant line. A stationary point of a differentiable concave function is a global maximum: concavity gives \(\ell(\gamma)\leq\ell(\beta)+U(\beta)^\mathsf{T}(\gamma-\beta)\), and the last term is zero at a stationary point. Strict concavity rules out two distinct maxima.

Worked Example: Evaluating Score and Curvature

Take three observations with covariates \(x_i^\mathsf{T}=(1,0),(1,1),(1,2)\) and outcomes \(0,1,1\). Evaluate the score and Hessian at \(\beta_0=0\), \(\beta_1=\log 2\). The linear predictors are \(0,\log 2,2\log 2\), so the fitted probabilities are \(1/2,2/3,4/5\). The residuals are \(-1/2,1/3,1/5\).

The intercept score is \(-1/2+1/3+1/5=1/30\). The slope score is \(0(-1/2)+1(1/3)+2(1/5)=11/15\). The score is not zero, so these coefficients do not solve the likelihood equations.

The three curvature weights are \(1/4,2/9,4/25\). Therefore,

H(\beta)=- \begin{pmatrix} \frac{1}{4}+\frac{2}{9}+\frac{4}{25} & \frac{2}{9}+\frac{8}{25}\\[4pt] \frac{2}{9}+\frac{8}{25} & \frac{2}{9}+\frac{16}{25} \end{pmatrix} =- \begin{pmatrix} \frac{569}{900}&\frac{122}{225}\\[2pt] \frac{122}{225}&\frac{194}{225} \end{pmatrix}.

The negative sign and positive curvature weights exhibit the expected downward curvature. Because the intercept and predictor columns are linearly independent, the design matrix has full column rank and this Hessian is negative definite.

Common Mistakes

Keeping an extra variance factor in the score. A tempting but incorrect step is to multiply \(y_i-\pi_i\) by \(\pi_i(1-\pi_i)\) when differentiating the log-likelihood with respect to \(\beta\). That factor appears when differentiating \(\pi_i\) with respect to \(\eta_i\), but the derivative of the Bernoulli log-likelihood with respect to \(\eta_i\) simplifies to \(y_i-\pi_i\). The variance factor belongs in the Hessian.
Claiming that every residual must vanish. The score requires \(X^\mathsf{T}(y-\pi)=0\), not \(y_i=\pi_i\) for every observation. With binary outcomes, individual residuals are generally nonzero.
Assuming residuals sum to zero without an intercept. The unweighted sum is the score component for an intercept. If the model has no intercept, that balance is not generally imposed.
Treating a score root as guaranteed to exist. The derivative formulas hold for finite coefficients, but they do not alone prove that the likelihood has a finite maximizer or that the score equations have a finite root.

Key takeaway: For individual Bernoulli outcomes under a logit link, each score contribution is the covariate vector multiplied by the observed-minus-fitted residual. The Hessian is the negative weighted cross-product of the covariates, with weights \(\pi_i(1-\pi_i)\).

Check Your Understanding

Use the individual Bernoulli logit likelihood and the score and Hessian formulas to answer these questions.

Write the log-likelihood contribution from one observation with outcome zero and linear predictor \(\eta_i\).
For a model with an intercept, what does the intercept score equation say about the fitted residuals?
Why does \(\pi_i(1-\pi_i)\) appear in the Hessian but not as an extra multiplier in the score?
What property of the design matrix ensures that the Hessian is negative definite at finite coefficient values?
In a model without an intercept, is a zero sum of residuals required by the coefficient score equations? Explain.