Two algorithms for fitting constrained marginal models

We study in detail the two main algorithms which have been considered for fitting constrained marginal models to discrete data, one based on Lagrange multipliers and the other on a regression model. We show that the updates produced by the two methods are identical, but that the Lagrangian method is more efficient in the case of identically distributed observations. We provide a generalization of the regression algorithm for modelling the effect of exogenous individual-level covariates, a context in which the use of the Lagrangian algorithm would be infeasible for even moderate sample sizes. An extension of the method to likelihood-based estimation under $L_1$-penalties is also considered.

💡 Research Summary

This paper provides a thorough comparison of two principal algorithms for fitting constrained marginal models to discrete data: a Lagrange‑multiplier based method derived from Aitchison and Silvey (1958) and a regression‑based algorithm introduced by Colombi and Forcina (2001). Both approaches aim to maximize the multinomial log‑likelihood subject to a set of nonlinear constraints that are expressed in terms of marginal log‑linear parameters (MLLPs). The authors first review the definition of MLLPs, emphasizing the concepts of completeness and hierarchy, which guarantee that the model belongs to a curved exponential family and that the constraints are smooth.

The Lagrange‑multiplier algorithm formulates the constrained optimisation problem as the stationary points of the Lagrangian L(θ,λ)=l(θ)+λᵀh(θ), where h(θ)=Kᵀη(θ) represents the linear constraints on the marginal parameters. By expanding the score vector s(θ) and the constraint function h(θ) to first order around a current estimate θ₀ and substituting the expected information matrix F for the Hessian of the log‑likelihood, the authors derive the update equation (2). This equation is essentially a Newton step corrected by a term involving the Jacobian H of the constraints. The update can be written without explicitly solving for the Lagrange multipliers, which simplifies implementation. The authors note that step‑size adjustments may be required for stability, especially when the constraints are near the boundary of the parameter space.

The regression algorithm proceeds by constructing a design matrix X that spans the orthogonal complement of the constraint matrix K (so that KᵀX=0) and re‑parameterising the marginal parameters as η = Xβ. The transformation matrix R = ∂θ/∂ηᵀ is used to map the score and information from the canonical θ‑space to the η‑space, yielding transformed quantities (\bar s = Rᵀ s) and (\bar F = Rᵀ F R). A quadratic approximation of the log‑likelihood around θ₀ is then expressed as a function of β, leading to a weighted least‑squares problem. Solving this problem gives the β‑update (3), and back‑substituting through R provides the corresponding θ‑update (4). The authors prove that the combined steps (3)–(4) are mathematically equivalent to the single update (2) from the Lagrange‑multiplier method, establishing the equivalence of the two algorithms.

A detailed computational complexity analysis shows that the Lagrange‑multiplier method’s most demanding operation is the product (Kᵀ C , \text{diag}(Mπ)^{-1} M), which scales as O(r·u·t) where r is the number of constraints, u the dimension of the marginal parameter vector, and t the total number of cells. In contrast, the regression method requires computing the matrix R, an O(u·t² + t³) operation, making it less efficient in pure computational terms. However, the regression framework can be naturally extended to incorporate individual‑level covariates. By stacking subject‑specific design matrices (X_i) and treating each subject’s marginal parameters as η_i = X_i β, the regression algorithm avoids the explosive dimensionality that would arise if the Lagrange‑multiplier approach were applied directly (which would require a design matrix of size n·t for n subjects). This makes the regression method far more practical for moderate to large sample sizes with covariates.

The paper also discusses an extension to more general constraints of the form (h(θ)=A \log(Mπ)=0) where A is an arbitrary full‑row‑rank matrix. Although the Jacobian in this case does not simplify as in the homogeneous case, the same quadratic‑approximation strategy yields an update analogous to (3).

Further, the authors introduce an L₁‑penalised version of the regression algorithm. By adding a term λ∑|β_j| to the objective, the update becomes a coordinate‑wise soft‑thresholding step, enabling simultaneous variable selection and constraint satisfaction. This extension is particularly useful for high‑dimensional settings where sparsity is desired.

Convergence properties are examined: if the constraints are smooth (full‑rank Jacobian everywhere) and the algorithm converges to a point θ*, then θ* satisfies the constraints and is a stationary point of the constrained likelihood. The Karush‑Kuhn‑Tucker conditions guarantee that such a stationary point is a local maximum provided the observed information matrix with respect to β is positive definite. The authors advise checking eigenvalues of the observed information and using multiple starting values to guard against convergence to sub‑optimal local maxima, especially because the constrained likelihood may be non‑convex.

In summary, the paper demonstrates that both algorithms generate identical updates but differ in computational efficiency and scalability. The Lagrange‑multiplier method is faster for models with identical observations and modest numbers of constraints, whereas the regression algorithm is indispensable when individual‑level covariates or L₁ penalties are required. The work offers clear guidance on algorithm selection based on model complexity, sample size, and the need for sparsity, thereby providing a valuable reference for practitioners working with constrained marginal models.