Consistent and Scalable Composite Likelihood Estimation of Probit Models with Crossed Random Effects

Estimation of crossed random effects models commonly require computational costs that grow faster than linearly in the sample size $N$, often as fast as $\Omega(N^{3/2})$, making them unsuitable for large data sets. For non-Gaussian responses, integrating out the random effects to get a marginal likelihood brings significant challenges, especially for high dimensional integrals where the Laplace approximation might not be accurate. A formula that is consistent under exact integration may fail to yield consistent estimates when numerical integrations are used. We develop a composite likelihood approach to probit models that replaces the crossed random effects model by some hierarchical models that require only one dimensional integrals. We show how to consistently estimate the crossed effects model parameters from the hierarchical model fits, using recent developments in adaptive Gauss-Hermite quadrature. We prove that the computation scales linearly in the sample size. We illustrate the method on about five million observations from Stitch Fix where the crossed effects formulation would require an integral of dimension larger than $700{,}000$.