Estimation of crossed random effects models commonly incurs computational costs that grow faster than linearly in the sample size, often as fast as, making them unsuitable for large datasets. For non-Gaussian responses, integrating out the random effects to obtain a marginal likelihood poses significant challenges, especially for high-dimensional integrals for which the Laplace approximation may not be accurate. In this article 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. The method is illustrated by applying it to approximately five million observations from Stitch Fix, where the crossed effects formulation would require an integral of dimension larger than.
