ABSTRACT

Digital learning platforms enable large-scale randomized controlled trials (RCTs) in real classroom settings, where students are nested within classrooms and outcomes are shaped by both individual characteristics and shared environments. Although randomization at the student level helps balance these factors, differences across classrooms may still persist, motivating the need for covariate adjustment.

In this paper, we study whether and how accounting for classroom structure improves treatment effect estimation. Standard approaches, such as fixed- and random-effects models, adjust for classroom variation and can improve precision, but they rely on parametric assumptions and may fail to capture complex or nonlinear relationships.

We build on the Leave-One-Out Potential outcomes (LOOP) estimator, which uses machine learning to flexibly model outcomes while preserving unbiasedness. However, LOOP struggles with high-cardinality categorical variables, such as classroom identifiers. To address this limitation, we propose Ordered LOOP, a method designed to better incorporate such variables.

Through simulations and real data applications, we show that accounting for classroom structure improves precision. Ordered LOOP performs particularly well when outcome heterogeneity is complex or nonlinear, while remaining competitive in simpler settings. Overall, our results highlight the importance of modeling classroom structure and demonstrate that Ordered LOOP provides a flexible, robust, and unbiased approach for analyzing hierarchical educational data.

Keywords

1. INTRODUCTION

Digital learning platforms allow researchers to conduct large-scale A/B tests and collect rich student-level data in real classroom settings. As their use grows, so does the need for methods that deliver valid and precise causal estimates in real-world environments. A key challenge is the hierarchical structure of educational data: students are nested within classrooms, teachers, and schools, and outcomes reflect both individual characteristics and shared learning environments [2].

Although student-level randomization is intended to balance these factors, real-world experiments often fall short. For example, analyses of ASSISTments experiments [6] show that many classrooms include only a single participating student [11], limiting the effectiveness of randomization and motivating covariate adjustment. Prior work [2] shows that accounting for classroom structure using fixed- or random-effects models can improve precision. However, these approaches rely on parametric assumptions and may fail to capture complex, nonlinear, or interactive relationships, especially when group structures are unknown.

To address these limitations, we study the Leave-One-Out Potential outcomes (LOOP) estimator [4, 14], which uses machine learning to flexibly model outcomes while preserving unbiasedness. Random forests are well suited to this framework because they naturally produce out-of-bag predictions aligned with the leave-one-out structure. However, standard implementations struggle with high-cardinality categorical variables, such as classroom identifiers, due to computational constraints.

We propose Ordered LOOP, an extension that enables effective adjustment for such variables. We evaluate the method using simulations based on real-world data and real experimental datasets [12] . Our results show that accounting for classroom structure improves precision, and that Ordered LOOP performs particularly well when classroom effects are complex or nonlinear, while remaining competitive in simpler settings. In these complex settings, it outperforms common alternatives, including OLS, OLS with classroom fixed effects, and hierarchical linear models with random effects.

The remainder of the paper is organized as follows. The next section reviews the potential outcomes framework, LOOP, and Ordered LOOP. Section 3 describes the simulation design, Section 4 presents the results, and Section 5 concludes.

2. CAUSAL ESTIMATOR

Consider a randomized controlled trial with \(N\) units indexed by \(i = 1, \dots , N\). Treatment assignments \(tr_i \in \{0,1\}\) are independent across units, with \(P(tr_i = 1) = p_i \in (0,1)\). For each unit, we observe outcome \(Y_i\) and covariates \(Z_i\). Under the potential outcomes framework [8], unit \(i\) has potential outcomes \(t_i\) and \(c_i\), and the observed outcome is \(Y_i = tr_i t_i + (1 - tr_i)c_i\). The individual treatment effect is \(\tau _i = t_i - c_i\), and the average treatment effect is \(\bar {\tau } = \frac {1}{N} \sum _{i=1}^N \tau _i\).

2.1 Leave-One-Out Potential outcome

Define the counterfactual expected outcome as \(m_i = (1 - p_i)t_i + p_i c_i\), with estimate \(\hat {m}_i\). Let the inverse probability weight be \(U_i = \frac {1}{p_i}\) if \(tr_i = 1\) and \(U_i = \frac {1}{1 - p_i}\) if \(tr_i = 0\). The unit-level estimator \(\hat {\tau }_i = (Y_i - \hat {m}_i) U_i\) is unbiased for \(\tau _i\) when \(\hat {m}_i\) is independent of \(tr_i\) [14]. The average treatment effect is estimated by \(\hat {\tau } = \frac {1}{N} \sum _{i=1}^N \hat {\tau }_i\), and the estimator has bounded variance when \(\hat {m}_i \approx m_i\) [14].

In practice, \(\hat {m}_i\) is obtained via leave-one-out prediction: for each unit \(i\), a model trained on the remaining \(N-1\) units produces predictions \(\hat {t}_i\) and \(\hat {c}_i\), with \(\hat {m}_i = (1 - p_i)\hat {t}_i + p_i \hat {c}_i\). Flexible methods such as random forests can capture nonlinearities and interactions while preserving unbiasedness, and naturally provide out-of-bag predictions aligned with this procedure.

2.2 Ordered LOOP

Standard random forest implementations face significant computational challenges when handling high-cardinality categorical predictors. For a predictor with \(k\) categories, the classical splitting rule considers \(2^{k-1}-1\) possible partitions, which quickly becomes infeasible even for moderate \(k\). Unlike continuous or ordinal variables, these splits cannot be represented by a single threshold. In practice, implementations often encode category assignments using bit representations, imposing hard limits (e.g., 32 or 64 levels) on the number of supported categories.

A common solution is to impose an ordering on categorical levels. For regression and binary classification, ordering categories by their mean response or class proportion yields the same optimal split as exhaustive search, while reducing the number of candidate splits from \(2^{k-1}-1\) to \(k-1\) [1, 3, 7]. When this ordering is computed once and reused throughout tree construction, splits can be implemented as simple thresholds, leading to substantial computational gains with minimal loss in predictive performance [13].

We adopt the ranger implementation, which supports this approach via respect.unordered.factors = "order". In this setting, factor levels are ordered once based on their mean response (for regression) or outcome class proportions (for classification), and are not reordered at each split [5], where “class” refers to the response category rather than classroom groupings.

Incorporating ordered factors into LOOP introduces two challenges. First, the imputed regression function \(\hat {m}_i\) must remain independent of each observation’s outcome. To ensure this, we use \(K\)-fold cross-fitting: the data are randomly partitioned into \(K\) folds, and for each fold, separate models are trained on the remaining \(K-1\) folds for the treated and control groups. Out-of-fold predictions \(\widehat {\mu }_1(tr)=\mathbb {E}[Y(1)\mid tr]\) and \(\widehat {\mu }_0(tr)=\mathbb {E}[Y(0)\mid tr]\) are then obtained for the held-out observations.

Second, arm-specific \(K\)-fold training can produce categorical levels in the test fold that are absent from the corresponding training data, leading to unstable predictions under ordered encoding. To address this, we adopt a row-level fallback strategy: for observations with unseen levels, we use predictions from a numeric-only imputation model (standard LOOP), while retaining ordered-ranger predictions for all other observations.

This hybrid approach leverages high-cardinality categorical information while ensuring stability and preserving the flexibility and unbiasedness of the LOOP estimator. When high-cardinality categorical variables are informative about the outcome, Ordered LOOP can exploit the induced ordering structure to reduce variance, effectively lowering the standard error and improving precision relative to standard LOOP. In settings where such structure is weak or absent, its performance remains comparable to LOOP.

3. SIMULATION DESIGN

Evaluating treatment effect estimators requires ground truth, which is unavailable in real experiments. We therefore simulate realistic settings using remnant [4, 9, 10] logs, which consist of observational student–assignment records aligned with experimental contexts.

3.1 Focal Measures

We use three variables: (i) student performance (inverse mastery speed, bounded in \([0,1/3]\), with 0 indicating non-completion [2]), (ii) prior performance (prior percent correct), and (iii) class membership (de-identified classroom ID).

3.2 Data Generation

We generate synthetic A/B tests using a bootstrap-style resampling procedure applied to the remnant data. This approach preserves the empirical distributions of class sizes, class-level heterogeneity, and student-level noise, while allowing controlled specification of treatment effects and outcome models.

We consider two data-generating processes. First, we replicate the baseline design of Closser et al. [2]. Second, we construct a more realistic setting in which class membership interacts with covariates and induces nonlinear outcome relationships. In this setting, standard fixed- and random-effects models are misspecified, allowing us to assess whether flexible machine learning methods can better capture the underlying structure.

3.2.1 Linear relationship

Under the linear specification, outcomes are generated as

Parameters are calibrated from the remnant data. We fit three models relating inverse mastery speed to prior percent correct: (i) OLS, (ii) OLS with class fixed effects, and (iii) a mixed-effects model with class-level random intercepts. The class-specific intercepts from the fixed-effects model define \(\text {class\_coef}\), capturing empirical classroom heterogeneity without imposing distributional and independence assumptions. The coefficient on prior performance is set to the average estimate across the three models (0.16).

Student-level noise \(\eta \) is generated by resampling residuals from the mixed-effects model, which provides a parsimonious decomposition of within- and between-class variation.

Each simulated dataset is constructed by sampling \(n_{\text {classes}}\) classrooms with replacement from the remnant data and retaining all associated observations. For each classroom, we draw a class-specific effect by resampling from the empirical distribution of the estimated class-level coefficients \(\text {class\_coef}\) obtained from the observed data.

At the observation level, we resample \(\eta \) and assign treatment independently with probability 0.5. Outcomes are then generated according to Equation (1).

3.2.2 Nonlinear relationship

Under the nonlinear specification, we introduce a class-relative threshold effect,

This component captures class-relative dynamics in which outcomes depend on a student’s performance relative to the class baseline. Higher-performing students may experience smaller gains (e.g., ceiling effects), while lower-performing students may benefit more from instruction, inducing a nonlinear pull toward the class mean.

Overall, this simulation framework preserves key empirical features of the remnant data while enabling controlled comparisons of estimator performance under both linear and nonlinear outcome-generating mechanisms.

4. RESULTS

The models evaluated in this study are summarized in Table 1 and are selected to address our central question: whether and how treatment effect estimation should account for the clustering of students within classrooms. We consider three classes of methods. First, we include approaches that ignore classroom structure, namely the two-sample \(t\) test and standard LOOP; the latter cannot incorporate classroom indicators when the number of categories is large (e.g., more than \(53\) levels). Second, we examine parametric models with fixed or random intercepts. Third, we evaluate Ordered LOOP, a machine-learning approach that flexibly captures classroom heterogeneity. Our goal is to assess when Ordered LOOP improves upon these alternatives.

4.1 Simulation Studies

All simulations fix the true treatment effect at \(0.2\) and are repeated \(100\) times.

4.1.1 Linear setting.

Table 2 reports results under the linear data-generating process. All methods recover the true effect accurately, but differ in precision. Methods that account for classroom structure achieve smaller standard errors than those that ignore it.

The random-effects model achieves the lowest variance, as expected under correct model specification. Ordered LOOP ranks second, indicating that it remains competitive even when the true outcome model is linear and fixed- or random-effects models are correctly specified.

In contrast, the standard LOOP estimator with covariate adjustment performs worse than the baseline \(t\) test. To better understand this result, we examined the imputation step \(\hat m\) in LOOP and found that the random forest used for outcome prediction explains a negative proportion of the outcome variance (\(-25.02\%\)). This indicates that the prediction model fails to capture useful signal in this setting, which in turn explains why covariate adjustment through LOOP does not improve—and in fact degrades—performance relative to the simple \(t\) test.

Overall, in the linear scenario, our findings are consistent with those of Closser et al. [2]: methods that explicitly account for classroom structure outperform approaches that ignore group-level heterogeneity, even when treatment is randomized at the student level. At the same time, Ordered LOOP remains a competitive alternative without requiring parametric assumptions about classroom effects.

4.1.2 Nonlinear setting.

Table 3 summarizes results under the nonlinear data-generating process. In this setting, parametric models are misspecified, and performance differences become more pronounced.

All methods recover the true effect accurately again. Ordered LOOP achieves the lowest standard error, followed by standard LOOP, demonstrating the advantage of flexible machine learning methods in capturing nonlinear and interaction effects. In contrast, fixed- and random-effects models are less efficient due to model misspecification.

Notably, standard LOOP performs well despite not explicitly including classroom indicators, suggesting that flexible modeling can partially compensate for missing group structure. However, incorporating classroom information through Ordered LOOP yields additional gains.

Across both settings, accounting for classroom structure consistently improves precision. Ordered LOOP performs robustly, remaining competitive in linear settings and outperforming alternatives in nonlinear settings, highlighting its effectiveness across both simple and complex outcome-generating processes.

4.2 Real-World Applications

4.2.1 Remnant data application

We construct a semi-synthetic dataset by assigning treatment at random, \(tr_i \sim \mathrm {Bernoulli}(0.5)\), since the remnant data contain no real treatment. Results are reported in Table 4. As expected, all estimates are close to zero.

Table 4 shows that methods accounting for classroom structure outperform those that ignore it: the hierarchical linear model achieves the smallest standard error, followed by Ordered LOOP and OLS with classroom fixed effects. Ordered LOOP is unbiased by construction, whereas the random-effects model relies on distributional assumptions that may not hold in practice.

4.2.2 LLM-based experiment.

We further evaluate all methods on a real experiment involving GPT-4–generated explanations [12]. Results are reported in Table 5.

The point estimates differ across methods, and because this is a real A/B test, the true effect is unknown. Theoretically, the t-test, fixed-effects model, LOOP, and Ordered LOOP are unbiased, while the random-effects model may be biased under misspecification.

Accounting for classroom clustering generally improves precision, but model misspecification can substantially increase variance. Ordered LOOP achieves the lowest standard error, outperforming the \(t\)-test baseline (0.0078 vs. 0.0083). In contrast, the random-effects model exhibits the largest variance, with a standard error 44% higher than that of Ordered LOOP (0.0112 vs. 0.0078), likely due to misspecification of its underlying assumptions. This illustrates a key limitation of parametric approaches: when the true data-generating process is unknown, misspecification can reduce robustness and lead to poorer performance. In contrast, more flexible, machine-learning–based methods that rely less on parametric assumptions tend to produce more stable and efficient estimates in such settings.

Across simulations and real data, two patterns emerge: accounting for classroom structure improves efficiency, and Ordered LOOP provides a robust and flexible approach that matches correctly specified parametric models and outperforms them when such models are difficult to specify. We therefore recommend Ordered LOOP for reliable performance across a wide range of settings.

5. DISCUSSION AND CONCLUSION

We study treatment effect estimation in online educational experiments with student-level randomization, focusing on classroom-level heterogeneity. Across simulations and real-world applications, methods that account for classroom structure, including fixed effects, random effects, and Ordered LOOP, consistently outperform those that ignore it, indicating that randomization alone may be insufficient for precise estimation.

When classroom effects are approximately linear, fixed- and random-effects models perform best under correct specification [2]. However, when effects are nonlinear or interact with covariates, parametric models may be misspecified and lose efficiency. In such settings, Ordered LOOP flexibly captures complex relationships while remaining competitive in simpler cases and retaining interpretability for estimating the average treatment effect.

Overall, accounting for clustering is essential for precision. Ordered LOOP offers a robust and flexible approach without requiring strong parametric assumptions. Future work could combine LOOP-style imputation with parametric models to further improve efficiency and interpretability.

6. ACKNOWLEDGMENTS

The research reported here was supported by the Institute of Education Sciences, U.S. Department of Education, through Grant R305D210031. The opinions expressed are those of the authors and do not represent views of the Institute or the U.S. Department of Education.

7. REFERENCES

1. L. Breiman, J. Friedman, R. A. Olshen, and C. J. Stone. Classification and regression trees. Chapman and Hall/CRC, 2017.
2. A. H. Closser, A. Sales, and A. F. Botelho. Should we account for classrooms? analyzing online experimental data with student-level randomization. Educational technology research and development, 72(5):2865–2894, 2024.
3. W. D. Fisher. On grouping for maximum homogeneity. Journal of the American statistical Association, 53(284):789–798, 1958.
4. J. A. Gagnon-Bartsch, A. C. Sales, E. Wu, A. F. Botelho, J. A. Erickson, L. W. Miratrix, and N. T. Heffernan. Precise unbiased estimation in randomized experiments using auxiliary observational data. Journal of Causal Inference, 11(1):20220011, 2023.
5. T. Hastie, R. Tibshirani, J. Friedman, et al. The elements of statistical learning, 2009.
6. N. T. Heffernan and C. L. Heffernan. The assistments ecosystem: Building a platform that brings scientists and teachers together for minimally invasive research on human learning and teaching. International Journal of Artificial Intelligence in Education, 24(4):470–497, 2014.
7. B. D. Ripley. Pattern recognition and neural networks. Cambridge university press, 2007.
8. D. B. Rubin. Estimating causal effects of treatments in randomized and nonrandomized studies. Journal of educational Psychology, 66(5):688, 1974.
9. A. Sales, A. Botelho, T. Patikorn, and N. T. Heffernan. Using big data to sharpen design-based inference in a/b tests. In Proceedings of the eleventh international conference on educational data mining, 2018.
10. A. C. Sales, B. B. Hansen, and B. Rowan. Rebar: Reinforcing a matching estimator with predictions from high-dimensional covariates. Journal of Educational and Behavioral Statistics, 43(1):3–31, 2018.
11. D. Selent, T. Patikorn, and N. Heffernan. Assistments dataset from multiple randomized controlled experiments. In Proceedings of the Third (2016) ACM Conference on Learning@ Scale, pages 181–184, 2016.
12. E. Worden, K. Vanacore, A. Haim, and N. Heffernan. Scaling effective ai-generated explanations for middle school mathematics in online learning platforms. In Proceedings of the Twelfth ACM Conference on Learning@ Scale, pages 40–49, 2025.
13. M. N. Wright and I. R. König. Splitting on categorical predictors in random forests. PeerJ, 7:e6339, 2019.
14. E. Wu and J. A. Gagnon-Bartsch. The loop estimator: Adjusting for covariates in randomized experiments. Evaluation review, 42(4):458–488, 2018.

© 2026 Copyright is held by the author(s). This work is distributed under the Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International (CC BY-NC-ND 4.0) license.