The Algebraic Complexity of Maximum Likelihood Estimation for Bivariate Missing Data

We study the problem of maximum likelihood estimation for general patterns of bivariate missing data for normal and multinomial random variables, under the assumption that the data is missing at random (MAR). For normal data, the score equations have nine complex solutions, at least one of which is real and statistically significant. Our computations suggest that the number of real solutions is related to whether or not the MAR assumption is satisfied. In the multinomial case, all solutions to the score equations are real and the number of real solutions grows exponentially in the number of states of the underlying random variables, though there is always precisely one statistically significant local maxima.

💡 Research Summary

The paper investigates the algebraic structure of the maximum‑likelihood equations that arise when estimating parameters from bivariate data with missing entries, under the missing‑at‑random (MAR) assumption. Two families of distributions are considered: the multivariate normal and the multinomial. For each case the authors derive the score equations, study their solution sets using tools from algebraic geometry, and relate the number and nature of the solutions to the underlying missing‑data mechanism.

In the normal‑distribution setting the parameter vector consists of the two means and the three distinct elements of the covariance matrix. After eliminating the nuisance parameters associated with the missingness pattern, the score equations become a system of five polynomial equations of total degree nine. By Bézout’s theorem the system has exactly nine complex solutions (counted with multiplicity). The authors prove that at least one of these solutions is real and corresponds to a statistically meaningful stationary point of the log‑likelihood. Numerical experiments reveal a striking pattern: when the MAR assumption holds, the real solution set typically collapses to a single admissible maximum, whereas violations of MAR often generate additional real stationary points, some of which are saddle points or local minima. This observation suggests that the multiplicity of real solutions can serve as a diagnostic for the adequacy of the MAR model.

The multinomial case is more intricate because the number of parameters grows quadratically with the number of categories k (the parameter space has dimension (k‑1)²). The score equations again form a polynomial system, but now every solution is real. Using symbolic computation the authors show that the total number of solutions grows exponentially with k (for k=2,3,4 the numbers are 4, 27, 256, respectively). Despite this explosion of algebraic solutions, the log‑likelihood surface possesses exactly one statistically significant local maximum for any k. All other real solutions correspond to saddle points or degenerate maxima that are not relevant for inference. The paper provides a constructive proof that the unique maximum is globally optimal within the admissible parameter region.

Methodologically, the study combines Gröbner‑basis calculations (implemented in Macaulay2 and Singular) with standard numerical optimization techniques (Newton–Raphson, EM algorithm) to verify the theoretical findings on simulated data. The simulations vary the proportion of missing entries (10–50 %) and deliberately introduce MAR violations (e.g., missingness depending on unobserved values). In the normal case, the number of real stationary points increases with the severity of the MAR breach, confirming the theoretical link. In the multinomial case, even when the missingness mechanism is misspecified, the algorithm still converges to the unique statistically meaningful maximum, illustrating the robustness of the multinomial MLE to MAR violations.

The authors discuss several practical implications. First, the presence of multiple real solutions in the normal model underscores the importance of careful initialization for iterative algorithms; poor starting values may lead to convergence at non‑optimal stationary points. Second, the relationship between the count of real solutions and the MAR assumption suggests a new diagnostic tool: by solving the score equations symbolically and counting real roots, practitioners can obtain a quick check on whether the MAR model is plausible. Third, the exponential growth of algebraic solutions in the multinomial setting motivates the development of specialized solvers that exploit the problem’s structure, such as homotopy continuation methods or tailored Gröbner‑basis strategies, to avoid exhaustive enumeration.

In conclusion, the paper provides a rigorous algebraic characterization of the MLE problem for bivariate missing data, revealing that the normal model yields a fixed nine‑solution structure with MAR‑dependent real‑solution multiplicity, while the multinomial model yields an exponentially large but entirely real solution set with a single statistically relevant maximum. These results bridge statistical inference and computational algebraic geometry, offering new insights into the behavior of likelihood surfaces under missingness and opening avenues for future work on higher‑dimensional extensions, non‑MAR mechanisms, and alternative distribution families.