Scalable Estimation of Crossed Random Effects Models via Multi-way Discretization

Cross-classified data frequently arise in scientific fields such as education, healthcare, and social sciences. A common modeling strategy is to introduce crossed random effects within a regression framework. However, this approach often encounters serious computational bottlenecks, particularly for non-Gaussian outcomes. In this paper, we propose a scalable and flexible method that approximates the distribution of each random effect by a discrete distribution, effectively partitioning the random effects into a finite number of representative values. This approximation allows us to express the model as a multi-way discrete structure, which can be efficiently estimated using a simple and fast iterative algorithm. The proposed method accommodates a wide range of outcome models and remains applicable even in settings with more than two-way cross-classification. We theoretically establish the consistency and asymptotic normality of the estimator under general settings of classification levels. Through simulation studies and real data applications, we demonstrate the practical performance of the proposed method in logistic, Poisson, and ordered probit regression models involving cross-classified structures.