From a Logistic Likelihood to a Numerical Algorithm
A logistic propensity model specifies treatment probabilities through covariates. Even when its likelihood equations have no closed-form solution, they can be solved numerically by repeatedly updating the coefficient estimates. Iteratively reweighted least squares (IRLS) gives each update as a weighted least-squares calculation. Understanding that calculation makes a software fit easier to check: the fitted values, weights, and update must all agree with the model.
Let \(A_i\) equal 1 if person \(i\) receives treatment and 0 otherwise, and let \(X_i\) be a column vector of covariates that includes an intercept. The logistic propensity model is
Here \(\beta\) is the coefficient vector and \(e_i(\beta)\) is the modeled probability of treatment for person \(i\). The model assumes that, conditional on the covariates, the treatment indicators are independent Bernoulli variables with these probabilities. The likelihood calculation also presumes the logistic mean model is the model being fitted. These are modeling assumptions; they do not, by themselves, establish causal identification or guarantee covariate balance.
For the usual unpenalized maximum-likelihood calculation, the design matrix must have full column rank, and the data must permit a finite maximum-likelihood estimate. If treatment status is perfectly or quasi-perfectly separated by the covariates, a finite estimate may not exist. IRLS can then produce ever-larger coefficients rather than converge to a finite solution.
The Score, Information, and Weighted Least Squares
For independent Bernoulli treatment indicators, the log likelihood, apart from constants, is
Differentiating with respect to \(\beta\) gives the score vector. Differentiating again gives the negative Hessian, called the information matrix in this logistic model. If \(X\) is the matrix with rows \(X_i^\mathsf{T}\), \(A\) is the vector of treatment indicators, and \(W\) is diagonal, then
In these expressions, \(e\) is the vector of fitted probabilities. A Newton or Fisher-scoring update solves \(I(\beta^{(t)})\Delta^{(t)}=U(\beta^{(t)})\), then sets \(\beta^{(t+1)}=\beta^{(t)}+\Delta^{(t)}\). In logistic regression the observed information and Fisher information have this same form.
The same update can be written as weighted least squares. At the current estimate, define the linear predictor \(\eta_i=X_i^\mathsf{T}\beta^{(t)}\) and the working response
The next estimate minimizes \(\sum_i W_{ii}(z_i-X_i^\mathsf{T}b)^2\) over \(b\). To see why, its normal equations are \(X^\mathsf{T}WXb=X^\mathsf{T}Wz\). Substituting the definition of \(z\) gives
Thus \(b=\beta^{(t)}+\Delta^{(t)}\), exactly the Newton update, provided the information matrix is invertible. The working response is a temporary quantity in the computation; it is not a newly observed treatment outcome. Its values and the weights change at every iteration.
Worked Example: One Intercept-Only Newton Update
Suppose a hypothetical dataset has five people, with two treated and three untreated. Fit an intercept-only logistic model, and suppose the current intercept is \(\beta^{(t)}=-0.4\). Using the rounded probability \(e=0.4013\), calculate one update.
The score for the intercept is the number treated minus the sum of the fitted probabilities: \(U=2-5(0.4013)=-0.0065\). The information is \(I=5(0.4013)(0.5987)=1.2013\). Therefore the increment is
The negative score says that at the current intercept the model slightly overestimates the observed treated proportion \(2/5=0.4\); the update decreases the intercept. The arithmetic here consistently uses the rounded probability stated above. Recomputing with unrounded values can change the last digits.
Worked Treatment-Assignment Iteration
Consider six hypothetical people, with one binary covariate \(x_i\). The data are summarized below; each row represents three people in that covariate group.
| Covariate \(x_i\) | Number treated \(A_i=1\) | Number untreated \(A_i=0\) |
|---|---|---|
| 0 | 2 | 1 |
| 1 | 1 | 2 |
Fit a logistic model with an intercept and \(x_i\), so that \(\eta_i=\beta_0+\beta_1x_i\). Start at \(\beta^{(0)}=(0,0)^\mathsf{T}\). The initial probability is \(0.5\) for every person, and every weight is \(0.5(1-0.5)=0.25\). The score is
There are three people in each covariate group, and the sum of \(x_i\) is 3. Consequently,
Solving \(I\Delta=U\) gives \(1.5\Delta_0+0.75\Delta_1=0\) and \(0.75\Delta_0+0.75\Delta_1=-0.5\). Subtracting the first equation from the second yields \(-0.75\Delta_0=-0.5\), so \(\Delta_0=2/3\), and then \(\Delta_1=-4/3\). Thus the first updated coefficients and predictors are
Applying the logistic function to these exact predictors gives fitted probabilities approximately \(0.6608\) for \(x=0\) and \(0.3392\) for \(x=1\). The treatment assignment is not perfectly predicted in either group, and both covariate groups contain treated and untreated people. The design matrix has full rank because the covariate varies and an intercept is included. These checks support carrying out the iteration; they do not yet show that the algorithm has converged.
Worked Example: Checking a Further Update
Continue the six-person example from \(\beta^{(1)}=(2/3,-4/3)^\mathsf{T}\). Let \(q=\operatorname{logit}^{-1}(2/3)\approx0.6608\) and \(r=\operatorname{logit}^{-1}(-2/3)\approx0.3392\). The score components are \(U_0=3-3q-3r=0\), since \(q+r=1\), and \(U_1=1-3r\approx-0.0177\).
Both groups have weight \(w=q(1-q)=r(1-r)\approx0.2242\), so the information matrix is \(w\begin{pmatrix}6&3\\3&3\end{pmatrix}\). Solving the two update equations gives \(\Delta_0=-U_1/(3w)\approx0.0264\) and \(\Delta_1=U_1/(1.5w)\approx-0.0527\). The next estimate is therefore approximately \((0.6930,-1.3861)^\mathsf{T}\). The increment is smaller than the first one, consistent with movement toward the finite maximum-likelihood estimate. A practical fit would continue checking the score and likelihood rather than stop solely because one update looks small.
Convergence Checks and Failure Modes
A sound numerical fit should be checked using more than a software message saying “converged.” Useful checks include the size of the score vector, the change in coefficients or log likelihood, and whether the information matrix can be solved stably. The fitted probabilities should also be inspected for values extremely close to 0 or 1. Those values can make weights \(e_i(1-e_i)\) very small, which may make updates unstable or signal separation.
Full column rank matters because otherwise different coefficient vectors can give the same linear predictors, and \(X^\mathsf{T}WX\) is singular. Even with full rank, nearly redundant columns or extremely different covariate scales can make the linear system poorly conditioned. A numerical routine may use step-halving or another safeguard if a full Newton step does not improve the likelihood. Such safeguards can stabilize computation, but they do not create a finite maximum if the data are separated.
- Check the model quantities. At each step, verify that probabilities lie between 0 and 1, weights equal \(e_i(1-e_i)\), and the working response uses the current—not previous—probabilities.
- Watch for separation. Coefficients that grow while probabilities approach 0 or 1 can indicate that no finite maximum exists.
- Distinguish optimization from model assessment. A small score indicates that the fitted likelihood equations are nearly solved. It does not prove the logistic specification is appropriate or establish balance on covariates.
Direct Newton updates and IRLS are algebraically equivalent in this model: one solves for a coefficient increment, while the other solves a weighted least-squares problem. General-purpose likelihood optimizers can instead use line searches or other numerical strategies and may be more robust in difficult cases. Whatever method is used, the target remains the same specified likelihood. Numerical convergence is not a substitute for considering model assumptions, and a fitted propensity score remains an estimate rather than a known treatment probability.
Check Your Understanding
Answer each question, then compare with the solutions below.
- In an intercept-only model with four treated and six untreated people, what is the score at \(\beta=0\)?
- Why does an IRLS working response change between iterations?
- For the six-person example, state the first updated linear predictor for each covariate group and its approximate fitted probability.
- What does a nearly zero score establish, and what does it not establish?
Solutions. (1) At \(\beta=0\), every fitted probability is \(0.5\), so \(U=4-10(0.5)=-1\). (2) The working response depends on \(\eta_i\), \(e_i\), and \(e_i(1-e_i)\), all of which are recalculated from the current coefficients. (3) The predictors are \(2/3\) for \(x=0\) and \(-2/3\) for \(x=1\); applying the logistic function gives approximately \(0.6608\) and \(0.3392\), respectively. (4) A nearly zero score indicates that the likelihood equations are nearly solved at the fitted coefficients. It does not establish that the model is correctly specified, that causal assumptions hold, or that covariates are balanced.