Computing score functions numerically using Mata

Specific econometric models - such as the Cox regression, conditional logistic regression, and panel-data models - have likelihood functions that do not meet the so-called linear-form requirement. That means that the model's overall log-likelihood function does not correspond to the sum of each observation's log-likelihood contribution. Stata's m1 command can fit said models using a particular group of evaluators: the d-family evaluators. Unfortunately, they have some limitations; one is that we cannot directly produce the score functions from the postestimation command predict. This missing feature triggers the need for tailored computational routines from developers that might need those functions to compute, for example, robust variance-covariance matrices. In this talk, I present a way to compute the score functions numerically using Mata's deriv() function with minimum extra programming other than the log-likelihood function. The procedure is exemplified by replicating the robust variance–covariance matrix produced by the clogit command using simulated data. The results show negligible numerical differences (e-09) between the clogit robust variance–covariance matrix and the numerically approximated one using Mata's deriv() function.