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

Estimation of crossed random effects models commonly requires 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. We develop a composite likelihood approach to probit models that replaces the crossed random effects model with 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. We find 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$.