Siyi Ruan
School of Statistics, East China Normal University, Shanghai, China.
DOI: 10.4236/ijis.2025.154012   PDF    HTML   XML   104 Downloads   611 Views   Citations

Abstract

Special education services are designed to provide tailored support for students with diverse learning needs, with the expectation of improving academic achievement. This study examines the causal effect of special education services on fifth-grade mathematics performance using data from the Early Childhood Longitudinal Study, Kindergarten Class of 1998-99 (ECLS-K). To address confounding, we employ both traditional and modern causal inference methods. First, propensity score matching (PSM) is implemented to create balanced treatment and control groups. Complementarily, three machine learning-based approaches—Bayesian Additive Regression Trees (BART), Causal Forests (CF), and Double Machine Learning (DML)—are applied to flexibly capture nonlinear relationships and treatment heterogeneity. Across all methods, we find consistently negative average treatment effects, with PSM and DML estimating stronger impacts, while BART and CF reveal moderately negative but significant effects. Importantly, machine learning models uncover substantial individual variation, indicating that while many students experience adverse outcomes, some benefit from targeted interventions. Variable importance analysis highlights early academic ability, socioeconomic background, and parental expectations as key predictors of later achievement, suggesting that observed disadvantages largely predate special education assignment. These findings underscore the methodological value of combining traditional and machine learning approaches for causal inference and carry important policy implications for designing equitable, preventive, and individualized educational interventions.

Keywords

Special Education, Causal Inference, Average Treatment Effect, Propensity Score Matching, Causal Machine Learning

Share and Cite:

1. Introduction

Special education services are the cornerstone of contemporary education policies, aimed at providing targeted assistance to students with different learning needs. Although implemented indiscriminately, their impact on student performance remains to be determined. Previous studies have produced conflicting evidence, often confused with selection bias, as students receiving special education often exhibit systematic differences in family background, cognitive abilities, and school environment compared to other students in the same school or student group. Therefore, the causal impact of special education is crucial in providing information for policies and ensuring effective and equitable intervention measures.

This study utilized data from the 1998-1999 Longitudinal Study of Preschool Classes (ECLS-K) to address the challenge of estimating the causal impact of special education services on fifth-grade mathematics grades. This dataset provides rich longitudinal information about academic performance, demographic background, family environment, and school environment, allowing for the control of all confounding factors.

To ensure control over selection bias, we have adopted both traditional and modern methods. Firstly, we adopt propensity score matching (PSM), a popular selection method in educational research, to create a balanced treatment group and an untreated group. Secondly, we supplement this with three machine learning supported causal inference methods—Bayesian Additive Regression Tree (BART), Causal Forest (CF), and Dual Machine Learning (DML)—which allow for flexible modeling of nonlinear relationships and heterogeneity in treatment outcomes. The similarity in this approach allows us to consider the distribution of average treatment outcomes (ATE) and individual treatment outcomes (ITEs), providing a more detailed perspective on the impact of special education.

Our empirical findings indicate consistently negative effects of special education on fifth-grade mathematics scores, with variation in magnitude across methods. Importantly, machine learning-based approaches reveal substantial heterogeneity, suggesting that while many students experience adverse outcomes, some may benefit from targeted interventions. Variable importance analysis further highlights early academic ability, socioeconomic status, and parental expectations as key predictors of later achievement, underscoring the role of pre-existing developmental and contextual factors in shaping educational outcomes.

Building on these considerations, this study pursues three specific objectives:

1) to clarify how different causal inference approaches shape our understanding of the effectiveness of special education.

2) to examine the extent to which treatment effects vary across individuals.

3) to identify the underlying factors that explain such heterogeneity.

In doing so, the analysis aims to provide both methodological insights for educational researchers and practical implications for the design of equitable and targeted special education policies.

This paper is organized as follows. Section 2 covers the literature and provides the theoretical background for our study. Section 3 describes the methodology, including both propensity score-based and machine learning-based estimation methods. Section 4 presents empirical results, ATE and ITE estimates, and variable importance analysis. Section 5 concludes with a discussion of the main results, their methodological and policy implications, and avenues for future research.

While multiple academic domains are available in the ECLS-K dataset, this study focuses on fifth-grade mathematics achievement since math is considered a major academic skill with strong predictive power for later educational attainment and labor market outcomes. Moreover, mathematics performance in fifth grade reflects cumulative learning from early schooling, making it an especially sensitive indicator of how special education services influence long-term academic trajectories.

2. Literature Review

Accurate estimation of treatment effects is a fundamental challenge in data science and causal inference. It is vital to understand how interventions impact outcomes to make good decisions in health and economics, among other fields.

2.1. Causal Inference Foundation

Causal inference (CI) is a core field of data science that is focused on inferring causal relationships through the use of counterfactual reasoning. CI provides a formal methodology for estimating and reasoning about the effect of interventions and presenting a systematic way of distinguishing causation from correlation [1] [2].

Crucial assumptions behind this framework are exchangeability [1], which states that treatment assignment is independent of potential outcomes given observed covariates:

$\left(\text{Y}\left(\text{1}\right)\text{,Y}\left(0\right)\right)$ $$ ⫫$$$\text{T}|\text{X}$

where Y(1) and Y(0) are the treatment and control potential outcomes, T is the indicator of treatment, and X is the vector of observed covariates. Exchangeability is intuitively the fact that treated and untreated groups become equivalent once we control for observed covariates as if they were randomized.

The Stable Unit Treatment Value Assumption (SUTVA) assumes that the potential outcomes of a unit are independent of other units’ treatment statuses [3] and that any particular treatment is in one, necessarily defined form. Mathematically,

$Y_i\left(t\right) \text{is well-defined and unaffected by }{T}_j\text{ for all }j\ne i$

These assumptions ensure that treatment effects are not confounded by external factors. Central to CI is the average treatment effect (ATE), which quantifies an intervention’s impact across a population:

$ATE=\mathbb{E}\left[Y\left(1\right)-Y\left(0\right)\right]$

This parameter is essential for evaluating the overall effectiveness of an intervention and for making evidence-based policy decisions [4].

2.2. Challenges in Causal Inference

While theoretically robust, CI faces some practical challenges. Causality cannot be observed but can only be inferred from correctly designed studies. It is inherently difficult to differentiate causation from correlation [2]. Standard statistical models, without a causal model such as the Neyman-Rubin model, are not sufficient for drawing conclusions regarding causality. A strict experimental or quasi-experimental design must be employed to eliminate confounding variables and receive correct treatment effect estimates. Moreover, traditional CI methods rely on parametric assumptions that may not hold when dealing with high-dimensional data having complex relations, and this leads to biased estimates [5] [6]. Assumptions like exchangeability and SUTVA may be violated in real-world situations due to interference across units or due to unobserved confounders.

2.3. Machine Learning and Causal Inference Integration

Machine learning (ML) offers powerful solutions to rectify some of the limitations of conventional CI procedures. Conventional regression models are bound to fail on complex, high-dimensional data and create biased estimates [7] [8]. ML, in contrast, offers concise predictive advantages, typically besting conventional statistical methods and offering insight beyond conventional models [9]-[11]. Under exogeneity, the performance of ML has been found to be superior, which is why it has been implemented in data science. In the estimation of treatment effects, ML techniques can potentially significantly reduce the bias [12]. For instance, Opper [13] introduced a novel ML technique that improves the accuracy of ATE estimation in small-scale randomized controlled trials, pointing out the ability of ML to increase estimation precision in scenarios with small sample sizes. Beyond improved accuracy, ML methods are particularly well-suited for detecting heterogeneous treatment effects. Unlike traditional regression or matching approaches, they can flexibly capture nonlinear relationships, automatically select important covariates, and adapt to high-dimensional data structures without strong parametric assumptions. These features make them especially valuable in educational settings, where student responses to interventions are rarely uniform.

2.4. Limitations of Machine Learning

Despite these successes, ML models have limitations that must be recognized in order to develop more robust and interpretable models. While ML works well in environments similar to where it was trained, it fails on novel or diverse datasets. This limitation is due to models overfitting to spurious correlations in training data, which undermines their ability to generalize as well as to pick up on latent causal mechanisms [14] [15]. Furthermore, the majority of ML models, and especially deep neural networks, are “black boxes”, meaning it is hard to understand the decision-making process. This uninterpretability is particularly troubling in areas where decision rationale is important. Without transparency, stakeholders cannot validate or critically assess model outputs, and trust in these systems is limited [16].

2.5. Emergence of Causal Machine Learning

Classical ML is focused on the identification of patterns and associations but falls short in identifying cause-and-effect relations [17]. This deficiency hinders its operation in those tasks that require close familiarity with the way interventions or changes might affect results. Addressing this gap involves the integration of methods that are more than pattern finding, with an effort to develop models capable of reasoning and understanding causal relationships.

Causal Machine Learning (CML) unites ML techniques with CI in order to take advantage of both approaches’ strengths. CML enhances the estimation of individualized treatment effects (ITE) and enables personalized prediction of the outcome under different interventions [18].

These classical methods normally involve parametric assumptions that are not necessarily found in practice. In order to sidestep this limitation, researchers have turned to more versatile methods of modeling. Bayesian Additive Regression Trees (BART) is a prime example because it provides a nonparametric Bayesian approach to modeling the causal inference response surface. BART has also been discovered to outperform traditional methods like propensity score matching and linear regression, particularly in nonlinear settings, producing more accurate estimates of average treatment effects [19].

2.6. Theoretical Perspectives on Special-Education Effects

Understanding why special-education interventions are at times harmful or unwanted is informed by broad-based educational theories:

The Matthew effect indicates that early skill sets, such as reading, compound on themselves—those who advance earlier are able to continue advancing, and students who struggle fall increasingly behind if problems still exist with the lowest skills [20]. For special-education students learning in special-education environments, early academic lags trap them in a “poor-get-poorer” feedback loop.

In contrast, compensatory education theory proposes that focused interventions among under-resourced youth will bridge performance gaps. However, if the assistance is poor or inadequately designed, intervention will bolster deficit narratives and neglect structural inequities, thus undermining educational equity.

Furthermore, stigmatization and labeling can potentially undermine student self-esteem and motivation: being labeled “special needs” may trigger self-fulfilling prophecies that weaken the psychological conditions for effective learning.

2.7. Applications of Causal Machine Learning

The fusion has given rise to real-world applications in various fields. In medicine, CML has been used to mimic complex, high-dimensional clinical data, allowing for more individualized and better treatment recommendations [21]. In economic policy evaluation, Belloni et al. [6] demonstrate how CML methods can effectively handle high-dimensional data to estimate treatment effects, as demonstrated in their study of the impact of 401 (k) plan participation on accumulation. Also, CML has been employed to optimize development aid allocation by forecasting heterogeneous treatment effects. A framework for optimizing aid distribution policies, particularly against HIV/AIDS, was suggested by Kuzmanovic et al. [22] using representation learning and counterfactual generation. Examples of applications illustrate the flexibility and power of CML in addressing actual issues.

In recent years, empirical research has begun to examine the heterogeneous impacts of special education services. Hwang [23] used Indiana administrative panel data and found that special education boosts mathematics and English achievement while increasing suspensions and absences, which indicates trade-offs between academic gains and behavioral costs. Hornby et al. [24] demonstrate that under specific instructional and support conditions, inclusive practices can yield measurable gains. Stark [25] calls for future research to deeply integrate individual-level heterogeneity, structural inequalities, and systemic supports in analyzing special education outcomes. Comparative work such as Dell’Anna et al. [26] further highlights that institutional context matters: while inclusive practices tend to show positive effects overall, the magnitude and direction vary across settings and frameworks.

Together, these theoretical and empirical insights underscore a multifactorial logic: special-education interventions do not exert uniform effects across all students or contexts. This is where methods like CML (causal machine learning) align well: beyond estimating an average treatment effect, CML helps identify subgroups or conditions.

3. Method

This section outlines the generalized workflow for estimating treatment effects from observational data, encompassing two primary analytical approaches: the Traditional Model-Based Approach and the Machine Learning-Based Approach.

3.1. Rubin Causal Model (RCM)

The causal effect estimation in this research relies on the Rubin Causal Model (RCM) method. Under the RCM method, each unit is capable of experiencing two possible outcomes: receiving treatment or not receiving treatment. There exists a single possible observable outcome, and thus it is not possible to observe the individual treatment effects. Therefore, the focus is on estimating the Average Treatment Effect (ATE). It calculates the difference between the expected outcome under the two situations.

Let $T\in \left\{0,1\right\}$ denote the binary treatment indicator, where $T=1$ indicates receipt of special education services and $T=0$ otherwise. $\text{Y}\left(1\right)$ and $\text{Y}\left(0\right)$ represent the potential outcomes under treatment and control, respectively. ATE is then defined as:

$\tau =\mathbb{E}\left[\text{Y}\left(1\right)-\text{Y}\left(0\right)\right]$

where:

$Y\left(1\right)$ : Potential outcome if the individual receives the treatment.

$Y\left(0\right)$ : Potential outcome if the individual does not receive the treatment.

$\mathbb{E}\left[\text{Y}\left(1\right)-\text{Y}\left(0\right)\right]$ : Expectation operator, representing the average over the population.

3.2. Baseline Method: Propensity Score Matching

In order to establish a baseline for causal estimation, we utilized Propensity Score Matching (PSM). PSM attempts to replicate the conditions of a randomized controlled trial by balancing the treated and untreated on observed characteristics. The method attempts to replicate a randomized experimental design by matching treated and untreated students with regard to similar observed characteristics, captured through their estimated probability of treatment [27].

We first estimated each student’s propensity score $\pi \left(X_i\right)=P\left(T_i=1 |X_i\right)$ using a logistic regression model:

$logit\left(\pi \left(X_i\right)\right)={\alpha }_0+X_i{}^{\top }\alpha$

where $T\in \left\{0,1\right\}$ denotes treatment assignment and $X_i$ is a vector of observed covariates.

Based on these scores, we conducted 1:1 nearest neighbor matching without replacement, using the logit-transformed score as the distance metric. This yielded a matched sample where treated and control groups are balanced on observed covariates.

To estimate the ATE, we ran a linear regression on the matched sample:

$Y_i=\alpha +\tau T_i+{ϵ}_i$

Then ATE was calculated as:

$ATE=\widehat{\tau }$

where $\widehat{\tau }$ is the estimated coefficient on the treatment variable $T_i$ . We computed heteroskedasticity-robust standard errors to ensure valid inference.

While PSM provides a transparent and interpretable framework for causal estimation, its reliance on parametric modeling of the treatment assignment mechanism poses potential limitations in high-dimensional or nonlinear settings. Logistic regression may fail to fully capture the complexity of real-world treatment selection. Therefore, we complement traditional approaches with modern machine learning-based estimators that flexibly account for confounding while preserving causal identification guarantees.

3.3. ML-Based Causal Inference Methods

3.3.1. Bayesian Additive Regression Trees

Bayesian Additive Regression Trees (BART) is a nonparametric Bayesian approach that models the outcome variable as a sum of regression trees, each capturing distinct patterns in the data. Intuitively, BART can be viewed as averaging over many small decision trees, each capturing a simple pattern in the data. By combining these trees within a Bayesian framework, BART produces flexible predictions and naturally accounts for uncertainty, much like consulting many different experts and then pooling their opinions with a measure of confidence.

This ensemble method naturally captures complex, nonlinear relationships and interactions without requiring manual specification, making it well-suited for high-dimensional settings where the true functional form is unknown [28]. Furthermore, BART’s Bayesian framework provides principled uncertainty quantification through posterior distributions, yielding credible intervals for treatment effect estimates, which is crucial for inference in observational studies [19]. BART inherently performs variable selection by determining the inclusion probabilities of covariates in the model, aiding in identifying relevant predictors without explicit feature selection procedures.

Within the potential outcomes framework, BART estimates the conditional expectation of the outcome given covariates and treatment:

$Y_i\left(t\right)={\displaystyle \sum }_{j=1}^m{g}_j\left(X_i,t\right)+{\epsilon }_i,{\epsilon }_i\sim N\left(0,{\text{σ}}^2\right),t\in \left\{0,1\right\}$

$\mathbb{E}\left[ Y|T,X\right]={\displaystyle \sum }_{t=1}^T{\displaystyle \sum }_{r=1}^R{f}_t\left(X;T\right)$

where:

$f_t\left(X;T\right)$ is the prediction from the $t^{th}$ regression tree for covariates $X$ and treatment $T$ .

The individual treatment effect (ITE) is computed as:

$ITE=\mathbb{E}\left[Y\left(1\right)|X\right]-\mathbb{E}\left[Y\left(0\right)|X\right]$

and the ATE is obtained by averaging ITEs:

$ATE=\frac{1}{N}{\displaystyle \sum }_{i=1}^{N}IT{E}_i$

In our implementation, BART models were fitted separately for treated and control groups. We initialized the model with a predefined number of trees and priors controlling tree depth and leaf parameter shrinkage. Posterior predictive distributions were generated for all units under both treatment and control conditions. Individual treatment effects were then calculated as the difference between predicted outcomes under treatment and control, and the ATE was computed as the mean of these ITEs. Posterior draws further allowed estimation of standard errors and 95% credible intervals, providing principled uncertainty quantification. In our implementation, the number of trees and prior shrinkage parameters followed standard defaults; sensitivity checks with different tree numbers produced consistent results, suggesting limited sensitivity to hyperparameter choice.

BART thus offers distinct advantages: its flexible nonparametric approach captures complex nonlinear relationships and interactions without manual specification, and its Bayesian framework naturally provides credible intervals for ATE estimates while allowing variable selection through inclusion probabilities.

3.3.2. Causal Forest (CF)

Causal Forest (CF) is a machine learning method designed specifically for causal inference in observational studies. Developed as an extension of Random Forests (RF), it adapts the principles of ensemble learning to estimate causal effects by leveraging recursive partitioning of the covariate space.

CF provides a flexible, nonparametric approach that can naturally handle complex, nonlinear relationships and interactions among covariates without requiring manual specification. This flexibility is especially beneficial in high-dimensional settings where the true functional form is unknown [29]. Secondly, CF employs an “honest” estimation strategy, wherein separate subsamples are used for tree construction and treatment effect estimation. This approach helps mitigate overfitting and ensures unbiased estimates [29].

The core idea of CF is to create a weighted similarity metric between observations, where similarity is defined not just by covariate values but by how those covariates influence the treatment effect. Here, it uses decision trees to maximize between-group heterogeneity in treatment effects in order to partition the data into subgroups.

$\text{Δ}\left(L\right)=\text{Var}\left({\left\{\widehat{{\text{τ}}_i}\right\}}_{i\in L}\right)-{\displaystyle \sum }_{k\in \left\{0,1\right\}}\frac{\left|L_k\right|}{\left|L\right|}\text{Var}\left({\left\{\widehat{{\text{τ}}_i}\right\}}_{i\in L_k}\right),$

where $L$ is a node and $L_k$ is its children.

Within each subgroup, it computes the difference in average outcomes between treated and untreated units via weighted averaging, akin to mini-randomized experiments.

$\widehat{\text{ITE}}\left(x\right)={\displaystyle \sum }_{i=1}^n{\text{α}}_i\left(x\right)\left(Y_i-\widehat{\text{μ}}\left(X_i\right)\right),{\text{α}}_i\left(x\right)=\frac{I\left(X_i\in \text{Leaf}\left(x\right)\right)}{\left|\text{Leaf}\left(x\right)\right|}.$

Then, the subgroup estimates are combined into a single ATE, which is weighted by the relevance of each subgroup to the overall population.

$ATE=\frac{1}{N}{\displaystyle \sum }_{i=1}^{N}IT{E}_i$

CF naturally identifies heterogeneous treatment effects across subgroups, while its built-in honesty mechanism, achieved through sample splitting, ensures more reliable and less biased effect estimates compared to conventional methods. The model was trained using default parameters (number of trees is 2000, minimum leaf size is 5), and robustness checks with alternative tree numbers and leaf sizes yielded qualitatively similar estimates, indicating stability of the results.

3.3.3. Double Machine Learning

Double Machine Learning (DML) is a semiparametric method. By modeling and removing covariate effects from both treatment and outcome variables, the causal relationship is isolated while high-dimensional confounders are handled. Conceptually, DML works by first removing (or “partialing out”) the parts of the outcome and treatment that can be explained by background covariates, and then estimating the treatment effect on the residual variation. This two-step process is akin to first leveling the playing field by adjusting for confounders, and then measuring the treatment’s unique contribution more cleanly.

DML utilizes Neyman orthogonality, which ensures that small estimation errors in the nuisance parameters have minimal impact on the estimation of the treatment effect, thereby reducing bias [30]. Additionally, it incorporates a cross-fitting strategy, where the data is split into folds, and nuisance parameters are estimated on separate folds from those used to estimate the treatment effect. It mitigates overfitting and enhances the robustness of the causal estimates.

DML creates orthogonalized versions of the treatment and outcome by removing their dependence on observed covariates, using cross-fitting to prevent overfitting in the nuisance parameter estimation. We respectively make the estimations based solely on background characteristics:

$\widehat{g}\left(X\right)=\mathbb{E}\left(Y|X\right)$

$\widehat{m}\left(X\right)=\mathbb{E}\left(T|X\right)$

The residual captures the unexplained variation:

$\tilde{Y}=Y-\widehat{g}\left(X\right)$

$\tilde{T}=T-\widehat{m}\left(X\right)$

with these purified variables, we estimate the ATE by:

$ATE=\frac{Cor\left({\tilde{T}}_i,{\tilde{Y}}_i\right)}{Var\left({\tilde{T}}_i\right)}=\frac{{{\displaystyle \sum }}_{i=1}^N{\tilde{T}}_i{\tilde{Y}}_i}{{{\displaystyle \sum }}_{i=1}^N{\tilde{T}}_i{}^2}$

This estimates how much unpredicted variation in math scores is associated with unpredicted variation in treatment assignment.

DML’s doubly robust property ensures reliable causal estimates even when only one of the models is correctly specified, while effectively handling the study’s high-dimensional nature by accommodating numerous covariates. We used random forests (500 trees, max depth = 5) to estimate nuisance functions, with 5-fold cross-fitting to mitigate overfitting; results were stable across variations in tree number and depth.

We employ BART, CF, and DML as complementary approaches to overcome the limitations of traditional methods. Together, they offer comprehensive causal insights—from overall effects to nuanced subgroup analyses.

3.4. Methodology

In this section, we describe the generalized workflow for estimating treatment effects from observational data, structured around two main analytical approaches: the Traditional Model-Based Approach and the Machine Learning-Based Approach. While the methodological assumptions have been detailed earlier, we focus on the computational procedures and data flow, starting from raw data input and moving to the final ATE estimation.

Let:

$X\in {ℝ}^{n\times p}$ : matrix of covariates

$T\in 0,1^n$ : binary treatment assignment vector

$Y\in {ℝ}^n$ : observed outcomes

Our objective is to estimate the ATE, denoted as:

$ATE=\frac{1}{n}{\displaystyle \sum }_{i=1}^n\left(Y_i\left(1\right)-Y_i\left(0\right)\right)$

where $Y_i\left(1\right)$ and $Y_i\left(0\right)$ are the potential outcomes under treatment and control, respectively.

3.4.1. Traditional Model-Based Approach

Step 1. Estimate Propensity Scores

Model the probability of treatment assignment given covariates using logistic regression:

$\widehat{e}\left(X_i\right)=P\left(T_i=1|X_i\right)$

Step 2. Perform Nearest Neighbor Matching

Using the estimated propensity scores, conduct 1:1 nearest neighbor matching without replacement. Matching is performed by minimizing the distance in logit-transformed scores:

$ℳ=\left\{\left(i,j\right):T_i=1,T_j=0,j=arg\underset{j\text{'}}{\min }\left|logit\left(\widehat{e}\left(X_i\right)\right)-logit\left(\widehat{e}\left(X_{j\text{'}}\right)\right)\right|\right\}$

where $ℳ$ denotes the matched sample of treated-control pairs.

Step 3. Estimate ATE

Compute the average difference in outcomes between matched treated and control units:

${\widehat{\tau }}_{PSM}=\frac{1}{ℳ}{\displaystyle \sum }_{\left(i,j\right)\in ℳ}\left(Y_i-Y_j\right)$

where $ℳ$ denotes the set of matched pairs, $m$ is the number of matched treated units, and $Y_i$ and $Y_j$ are the outcomes for the treated and matched control units, respectively.

3.4.2. Machine Learning-Based Approaches

To address limitations of traditional methods, we employ three machine learning techniques: Bayesian Additive Regression Trees (BART), Causal Forest (CF), and Double Machine Learning (DML).

1) Bayesian Additive Regression Trees (BART)

Step 1. Model Fitting

Fit BART models to predict outcomes under treatment and control:

$Y_i\left(t\right)={\displaystyle \sum }_{j=1}^m{g}_j\left(X_i,t\right)+{\epsilon }_i,{\epsilon }_i~N(0,{\sigma }^2),t\in \left\{0,1\right\}$

where $g_j$ represents the $j^{th}$ regression tree, and ${\epsilon }_i\sim N\left(0,{\text{σ}}^2\right)$ .

• Each tree captures distinct nonlinear patterns and interactions.

• Default hyperparameters are used (number of trees, tree depth, and prior parameters).

• A fixed random seed (2025) is set to ensure reproducibility.

• Posterior predictions are used to compute ITE and credible intervals.

Step 2. Estimate ITE

Compute the difference in predicted outcomes:

${\widehat{\tau }}_i={\widehat{Y}}_i\left(1\right)-{\widehat{Y}}_i\left(0\right)$

Step 3. Estimate ATE

Average the individual treatment effects:

${\widehat{\tau }}_{BART}=\frac{1}{n}{\displaystyle \sum }_{i=1}^n{\widehat{\tau }}_i$

2) Causal Forest (CF)

Step 1. Model Training

Train a causal forest to partition the covariate space and estimate treatment effects within each leaf node.

• Default hyperparameters are used (number of trees, minimum node size, sample fraction).

• Random seed is set to 2025 for reproducibility.

Step 2. Estimate ITE

For each individual, estimate the treatment effect based on their covariates:

${\widehat{\tau }}_i=\widehat{\tau }\left(X_i\right)$

Step 3. Estimate ATE

Average the estimated individual treatment effects:

${\widehat{\tau }}_{CF}=\frac{1}{n}{\displaystyle \sum }_{i=1}^n{\widehat{\tau }}_i$

3) Double Machine Learning (DML)

Step 1. Estimate Nuisance Functions

Use machine learning models to estimate the expected outcome and treatment assignment given covariates:

$\widehat{g}\left(X\right)=\mathbb{E}\left(Y |X \right),\widehat{m}\left(X\right)=\mathbb{E}\left(T|X \right)$

• Outcome model: Random Forest regression (num.trees = 500, max.depth = 5).

• Treatment model: Random Forest classification (num.trees = 500, max.depth = 5).

• 5-fold cross-fitting is applied.

• Random seed is set to 2025 to ensure reproducibility.

Step 2. Computer Residuals

Calculate residuals by removing the estimated effects of covariates:

$\tilde{Y}=Y-\widehat{g}\left(X\right),\tilde{T}=T-\widehat{m}\left(X\right)$

Step 3. Estimate ATE

Use the residuals to estimate the treatment effect:

${\widehat{\tau }}_{DML}=\frac{{{\displaystyle \sum }}_{i=1}^n{\tilde{T}}_i{\tilde{Y}}_i}{{{\displaystyle \sum }}_{i=1}^n{\tilde{T}}_i{}^2}$

3.5. Comparative Overview of Methods

To provide a high-level perspective on the methodological choices, Table 1 summarizes the main features, strengths, and limitations of the causal inference methods employed in this study. This overview allows readers to quickly grasp the methodological landscape before delving into the detailed descriptions of each approach.

Table 1. Comparative overview of causal inference methods.

Method	Model Type	Key Strengths	Limitations
PSM	Parametric, logistic regression	Transparent; interpretable; replicates randomized design	Sensitive to model misspecification; weak in high-dimensional or nonlinear settings
BART	Nonparametric, Bayesian trees	Captures nonlinearities; posterior uncertainty; implicit variable selection	Computationally intensive; tuning required; reduced interpretability
CF	Nonparametric, random forest extension	Handles high-dimensional data; detects heterogeneity; honesty reduces bias	Needs large samples; limited interpretability
DML	Semiparametric, orthogonalization	Doubly robust; high-dimensional confounders; cross-fitting reduces bias	Depends on ML models for nuisance; complex procedure

This overview facilitates a better understanding of the subsequent sections, where we present detailed implementation and estimation procedures for each method.

To further highlight the methodological distinction, Table 2 summarizes the analytical steps of the traditional propensity score-based approach compared with the machine learning-based methods. While PSM relies on explicitly modeling treatment assignment and matching individuals, the machine learning methods share a common workflow of model fitting, counterfactual prediction, and effect estimation, thus providing a more unified and flexible framework for causal inference.

Table 2. Step comparison between traditional and ML-based methods.

Approach	Main Steps
Traditional: PSM	1) Estimate propensity scores using logistic regression;
	2) Perform nearest-neighbor matching to construct balanced groups;
	3) Estimate ATE on the matched sample via regression.
Machine Learning: BART/CF/DML	1) Fit machine learning models to predict treatment and/or outcome;
	2) Generate counterfactual predictions or residualized variables;
	3) Derive individualized treatment effects (ITEs) and average treatment effect (ATE) estimates.

Figure 1. Unified workflow of the study. Both the traditional approach (PSM) and machine learning-based methods (BART, Causal Forest, and DML) follow a general three-step process: model estimation, counterfactual prediction, and treatment effect estimation.

All methods were applied to the same analytic sample after excluding cases with missing covariates. For propensity score matching (PSM), unmatched observations without suitable counterparts were discarded, which slightly reduced the sample size but maintained covariate balance. By contrast, the machine learning-based approaches retained the full analytic sample, with internal procedures such as sample splitting and honesty used to mitigate potential overlap issues. Consequently, while small differences in sample inclusion may contribute to variation in estimated effect magnitudes, the direction of effects remains robust across methods.

To assess robustness, we conducted supplementary checks varying model specifications and sample restrictions. For PSM, alternative matching algorithms (nearest neighbor vs. caliper matching) produced qualitatively similar estimates. For BART and CF, we varied the number of trees and prior/shrinkage parameters, while for DML, we adjusted the depth of the random forest learners. Across these variations, the direction and statistical significance of treatment effects remained consistent, and effect magnitudes differed only modestly. We also repeated analyses after trimming observations with extreme covariate values, which yielded comparable results.

To further aid reader comprehension, Figure 1 presents a unified workflow of the study.

4. Result

4.1. Data and Variable Description

This study utilizes data from the Early Childhood Longitudinal Study, Kindergarten Class of 1998-99 (ECLS-K), a nationally representative panel survey administered by the National Center for Education Statistics (NCES). The ECLS-K followed over 20,000 children from kindergarten entry in 1998-1999 through their elementary school years, collecting rich longitudinal information on their academic development, family environment, and school experiences.

All analyses were implemented in R. Specifically, we used the package MatchIt for propensity score matching, dbarts for Bayesian Additive Regression Trees (BART), grf for CF estimation, and DoubleML for Double Machine Learning. Data processing and visualization relied on the tidyverse (e.g., dplyr, ggplot2). Complete replication code is provided in the Appendix to ensure transparency and reproducibility.

For this analysis, we focus on observations from the spring of 2004, when most students were in the fifth grade. The treatment variable is F5SPECS, a binary indicator of whether a child received special education services (1 = Yes, 0 = No). The primary outcome is C6R4MSCL, the fifth-grade mathematics achievement score, scaled using Item Response Theory (IRT).

In this study, missing values in covariates were handled by a listwise deletion approach, which is commonly used in empirical practice. That is, any observation with missing values in one or more covariates used in the analysis was removed from the sample, yielding the final analytic sample of 7362 students. Among them, 429 students (5.8%) received special education services, while 6933 students (94.2%) formed the comparison group. Descriptive statistics of the outcome further revealed notable baseline differences: the control group had an average fifth-grade math score of 128.19 (SD = 22.67), whereas the treatment group averaged 108.97 (SD = 26.80). These baseline disparities highlight the need for rigorous causal inference methods to adjust for potential selection bias.

While the ECLS-K data provide cross-sectional and longitudinal sampling weights designed to recover nationally representative estimates, the current study does not apply these weights. This reflects our analytic focus: we aim to estimate conditional treatment effects within the analytic sample, rather than producing descriptive population means. Weighting is essential for unbiased descriptive statistics, but when covariates are explicitly included in the causal models, unweighted estimation can yield more efficient and unbiased treatment effect estimates [31] [32].

To address potential selection bias and confounding, a wide array of pre-treatment covariates is included in the analysis. These covariates were measured prior to the treatment and span the following domains:

• Demographic characteristics, including gender, race/ethnicity, and family socioeconomic status;

• Academic background, such as kindergarten-level reading and math achievement, and parental expectations for education;

• Family structure and early care environment, including single-parent status, parental age at first birth, number of siblings, and early childcare experiences;

• Health and developmental indicators, such as birth weight, fine and gross motor skills, and early intervention service usage;

• Behavioral and psychological traits, based on parent-rated attentiveness, impulsivity, verbal communication, and emotional functioning;

• Classroom and school-level aggregates, including average academic scores and contextual measures like classroom SES and racial composition.

Table 3. Summary statistics for selected analytical variables.

Variable	Description	Mean	SD	Min	Max
F5SPECS	Special Education Services (0 = No, 1 = Yes)	0.0583	0.2343	0	1
C6R4MSCL	Fifth Grade Math IRT Score	127.0713	23.3665	50.86	170.66
WKSESL	Socioeconomic Status (Continuous Index; Higher = Higher SES)	0.1666	0.7614	−4.75	2.75
MIRT	Kindergarten Math Score (IRT-Scaled)	33.3464	10.0866	11.85	98.99
P1EXPECT	Parental Educational Expectations (1 = ≤ High School, 6 = Graduate Degree)	4.1030	1.0336	1	6
P1DISABL	Parent-Reported Disability (0 = No, 1 = Yes)	0.1278	0.3339	0	1
ONEPARENT	One-Parent Household (0 = No, 1 = Yes)	0.1718	0.3773	0	1
wt_ounces	Birth Weight in Ounces	119.4999	20.5140	17	214
P1ATTENI	Attentiveness Rating (1 = Not Attentive, 4 = Highly Attentive)	1.8396	0.6202	1	4
Avg_MIRT	Class-Level Mean Math Score	33.3464	5.4657	16.0625	66.075

Table 3 below summarizes the descriptive statistics (mean, standard deviation, minimum, and maximum) for a selected set of key analytical variables.

To contextualize the analysis, Table 4 presents the same key variables disaggregated by treatment status (special education vs. non-special education), along with standardized differences. This comparison highlights baseline imbalances and motivates the need for causal adjustment methods.

Table 4. Baseline characteristics of treatment and control groups.

Variable	Treatment Mean (N = 429)	Control Mean (N = 6933)	Std. Diff.
Gender (1 = Female)	0.674	0.486	+0.38
White (1 = Yes)	0.818	0.746	+0.16
SES Index (WKSESL)	−0.041	0.179	−0.29
Kindergarten Math (MIRT)	26.72	33.76	−0.72
Kindergarten Reading (RIRT)	35.24	42.61	−0.65
Parent Expectation (P1EXPECT)	3.79	4.12	−0.31
One-Parent Household	0.219	0.169	+0.12
Birth Weight (oz)	117.34	119.63	−0.11
Disability (Parent Report)	0.382	0.112	+0.65
Grade 4 Math (C5R4MSCL)	85.75	103.93	−0.77
Grade 4 Reading (C5R4RSCL)	131.65	155.83	−0.84
Grade 5 Math (Outcome, C6R4MSCL)	108.97	128.19	−0.78

4.2. Propensity Score-Based Estimation

4.2.1. Propensity Score Estimation

To estimate the causal impact of special education services on fifth-grade mathematics achievement, we constructed propensity scores using logistic regression. The treatment indicator (F5SPECS) was regressed on a comprehensive set of pre-treatment covariates, including demographic characteristics, academic background, family structure, health indicators, behavioral traits, and classroom/school-level contextual measures.

The estimated propensity scores varied considerably across individuals, reflecting wide heterogeneity in the likelihood of receiving special education. As shown in Figure 2, treated students exhibited broader dispersion in their scores, while untreated students were heavily concentrated near zero. Despite this imbalance, the existence of common support between the two groups justified the application of PSM.

Common-support diagnostics revealed that treated students had propensity scores ranging from 0 to 0.928, while control students ranged from 0 to 0.768, yielding a common support interval of [0, 0.768]. Approximately 88% of the original sample was trimmed during matching, reflecting the limited overlap between treated and untreated groups. While this substantial reduction underscores the restrictiveness of PSM in this context, it also highlights the importance of complementing PSM with flexible machine learning approaches that retain the full sample.

Figure 2. Distribution of estimated propensity scores by treatment status prior to matching. While treated students show broader dispersion, both groups exhibit sufficient overlap, providing common support for PSM and justifying the matching procedure.

4.2.2. Covariate Balance Assessment

Figure 3. Covariate balance before and after PSM was measured using absolute standardized mean differences. Post-matching balance largely falls below 0.1, indicating that PSM effectively reduced observable differences and approximates a randomized setting.

We implemented one-to-one nearest neighbor matching without replacement, using the logit of the estimated propensity scores as the distance metric. To assess covariate balance before and after matching, we computed standardized mean differences (SMDs) for all pre-treatment covariates.

As displayed in Figure 3, the matching procedure successfully reduced pre-existing imbalances across most covariates. Post-matching SMDs were generally below 0.1, indicating adequate balance and lending support to the credibility of the matching design. This suggests that the matched sample approximates a randomized setting in terms of observed characteristics.

A full table of pre- and post-matching balance statistics for all covariates is provided in Appendix A, confirming that standardized mean differences fell below 0.1 across nearly all covariates.

4.2.3. Estimating Treatment Effects with PSM

Using the matched sample, we estimated the ATE of receiving special education services via linear regression. The outcome variable was the fifth-grade mathematics IRT score (C6R4MSCL), and the sole predictor was treatment status.

The estimated ATE was −6.819 (SE = 1.801, 95% CI: [−10.348, −3.289]), indicating a statistically significant and negative effect of receiving special education on math achievement. To account for the matched-pair structure, we also re-estimated the ATE using cluster-robust standard errors, yielding a nearly identical estimate (−6.819, SE = 1.649, 95% CI: [−10.051, −3.586]). This finding suggests that, conditional on observable characteristics, treated students performed worse in mathematics than their matched untreated peers.

This result is consistent with estimates from machine learning-based methods (see Section 4.3), reinforcing the interpretation that observed treatment assignments may reflect pre-existing disadvantage rather than effective intervention. In other words, special education services may not have fully offset the challenges that led to treatment assignments in the first place.

4.3. Machine Learning-Based Estimation

To validate and complement the PSM results, we implemented three machine learning methods—Bayesian Additive Regression Trees (BART), Causal Forest (CF), and Double Machine Learning (DML). These models offer flexible, data-adaptive approaches to causal estimation, especially in the presence of complex and high-dimensional confounding structures.

4.3.1. Estimating ATE with BART

The BART model was employed to estimate both the ATE and the individualized treatment effects (ITEs). It was trained using the same treatment indicator (F5SPECS), outcome variable (C6R4MSCL), and covariates as used in the previous model. Treated and untreated groups were modeled separately, allowing for imputation of counterfactual outcomes.

The estimated ATE was −4.115 (SE = 0.07, 95% CI: [−4.252, −3.977]), indicating a significant negative effect of special education services. The reported 95% intervals are Bayesian credible intervals derived directly from posterior draws of the BART model. The ITEs showed substantial heterogeneity, ranging from 30.250 to 17.645 (SD = 6.019), as depicted in Figure 4. While most students exhibited negative effects, some benefited mildly, highlighting variation in treatment response across individuals.

Figure 4. Distribution of estimated ITEs from BART. While most students experience negative impacts, the widespread evidence suggests substantial heterogeneity, with some students benefiting from special education.

4.3.2. Estimating ATE with CF

Figure 5. Distribution of ITEs from Causal Forest. Effects are predominantly negative and left-skewed, but with meaningful heterogeneity across individuals.

The Causal Forest (CF) model similarly used the full dataset and covariate set to estimate treatment effects. It yielded an ATE of −3.589 (SE = 0.865, 95% CI: [−5.284, −1.893]). Confidence intervals were obtained using the asymptotic variance estimates implemented in the grf package, which leverages honesty and sample-splitting to provide valid inference. The ITEs ranged from −14.299 to 0.240, with a standard deviation of 1.720. As shown in Figure 5, the ITE distribution was strongly left-skewed, reinforcing the presence of heterogeneous and predominantly negative effects.

The CF estimate aligns with BART in both direction and magnitude, and is consistent with the broader pattern of negative treatment effects observed across all methods. This discrepancy may reflect the machine learning models’ superior ability to capture nonlinearities and latent confounding, which standard regression-based methods overlook.

4.3.3. Estimating ATE with DML

The DML model used a partially linear regression framework, modeling both treatment and outcome flexibly with random forests. Standard errors were computed via the DoubleML package based on asymptotic theory, and 95% confidence intervals were constructed accordingly. The estimated ATE was −6.792 (SE = 1.068, 95% CI: [−8.886, −4.698]), the largest in magnitude among all methods. The ITEs were tightly concentrated (range: −3.560 to −0.052, SD = 0.502), as shown in Figure 6, suggesting a consistently negative treatment effect with limited individual variation.

This compact distribution likely stems from the PLR model’s structural assumptions, which constrain effect heterogeneity more than BART or CF.

Figure 6. Distribution of ITEs from DML. Effects are tightly clustered around a consistently negative value, reflecting limited heterogeneity under this specification.

4.3.4. Discussion of Heterogeneous Effects

Both BART and CF identified substantial heterogeneity in the impact of special education services. Socioeconomic status (SES), early mathematics readiness, and parental expectations emerged as the most influential moderators of treatment effects. Students from lower-SES backgrounds and those entering school with weaker foundational skills were disproportionately negatively affected, whereas relatively advantaged students showed smaller negative effects.

The convergence of results across BART and CF strengthens the credibility of these findings. Although DML does not provide variable importance measures in the same manner, subgroup analyses based on key covariates produced consistent patterns. Taken together, the evidence suggests that pre-existing disadvantages amplify the adverse impacts of special education services, highlighting the need for preventive and targeted interventions.

4.4. Summary of Results

Table 5. Summary of ATE and ITE across methods.

Method	ATE_Estimate	Std_Error	CI_Lower	CI_Upper	ITE_Min	ITE_Max	ITE_SD
PSM	−6.819	1.801	−10.348	−3.289
BART	−4.115	0.070	−4.252	−3.977	−30.250	17.645	6.019
CF	−3.589	0.865	−5.284	−1.893	−14.299	0.240	1.720
DML	−6.792	1.068	−8.886	−4.698	−3.560	−0.052	0.502

Table 5 presents a comparative summary of the ATE estimates and ITE variability across the four methods. As visualized in Figure 7, all four methods—PSM, BART, CF, and DML—consistently indicate a negative effect of receiving special education services on fifth-grade math performance. The PSM method estimated the most negative ATE (−6.819, SE = 1.801), closely aligned with DML (−6.792), while BART and CF yielded moderately negative effects (−4.054 and −3.589, respectively). These results suggest that, after adjusting for observed covariates, students who received special education services tended to underperform relative to their matched or modeled untreated counterparts.

However, the estimated magnitudes varied: PSM and DML produced larger negative effects, while BART and CF yielded more moderate but still significant effects. These discrepancies likely reflect methodological differences. PSM, while intuitive and transparent, depends on the successful balancing of observed covariates and may amplify residual imbalance, producing stronger negative estimates. DML, although designed to reduce bias through orthogonalization, relies on sample splitting and random forests for nuisance estimation, which can increase variability, but in our case, produced estimates close to PSM. By contrast, BART and Causal Forests leverage flexible, nonparametric models that capture nonlinearities and interactions, often leading to more conservative but robust estimates.

Taken together, these patterns underscore the trade-offs across methods: PSM and DML may exaggerate magnitudes if covariates remain imperfectly balanced, whereas BART and CF provide more stable estimates under weaker functional-form assumptions. The consistency in effect direction across all approaches strengthens confidence in the overall conclusion, while the methodological contrasts highlight the importance of interpreting effect sizes in light of each method’s assumptions and potential biases.

Figure 7. Estimated ATEs and 95% confidence intervals across PSM, BART, CF, and DML. All methods converge on negative effects, reinforcing robustness despite methodological differences.

4.5. Variable Importance and Contribution Analysis

To further enhance the interpretation of our findings and address the contribution of individual covariates, we conducted a variable importance analysis using a random forest model [33].

Figure 8. Variable importance analysis of math achievement predictors. Early academic ability, socioeconomic background, and parental expectations emerge as the strongest drivers of outcomes, highlighting pre-existing disadvantages among special education recipients.

The random forest model was trained on all pre-treatment covariates to predict students’ math scores. Variable importance was quantified using the percentage increase in mean squared error (%IncMSE) when each variable was permuted, a standard measure in tree-based models. Figure 8 displays the top 10 most important variables ranked by their predictive power.

The results highlight that early academic ability, particularly kindergarten math achievement (MIRT), was the strongest predictor of fifth-grade math performance. Socioeconomic status (WKSESL), parental educational expectations (P1EXPECT), and class-level academic environment (avg_MIRT) also emerged as major contributors. These findings reinforce the critical role of early cognitive development and family background in shaping long-term academic outcomes.

This analysis complements the causal estimates by identifying which baseline attributes most strongly account for variation in academic outcomes, offering useful insights for targeted early interventions.

5. Conclusions

This study compared propensity score matching (PSM) with machine learning approaches—BART, CF, and DML—in estimating the causal effect of special education on fifth-grade mathematics achievement. All methods revealed negative ATE, with PSM and DML showing stronger negative impacts, while BART and CF produced more moderate yet significant estimates. Machine learning methods additionally uncovered effect heterogeneity: although most students experienced adverse outcomes, certain subgroups appeared to benefit. Key moderators included early academic readiness, socioeconomic background, and parental expectations, indicating that disadvantages among special education students primarily stem from pre-existing factors.

These findings offer both methodological and practical implications. Methodologically, they demonstrate machine learning’s value in detecting heterogeneity that traditional methods might miss. Practically, they show that special education placement does not uniformly benefit students and may worsen disadvantages for those from low-SES or less-prepared backgrounds. This highlights the need for early interventions, individualized instruction, and preventive strategies such as expanded access to high-quality preschool and targeted academic support.

Several limitations warrant acknowledgment. First, the observational nature of data leaves potential for unobserved confounding. Second, measurement error in test scores and survey variables may attenuate effects. Third, focusing solely on fifth-grade mathematics limits generalizability to other subjects or grades. Fourth, each machine learning method presents specific challenges: BART’s computational intensity, CF’ increased variance from sample splitting, and DML’s dependence on proper nuisance model specification. More broadly, all methods rely on the untestable assumption of ignorability given observed covariates.

Future research could extend this work by incorporating longitudinal data to track dynamic effects across grades, examining subgroup heterogeneity by gender, race, or school context, and employing alternative causal inference strategies like instrumental variables or synthetic controls. Advances in deep learning estimators may also complement existing methods and offer new avenues for analyzing large-scale educational datasets.

Appendix A: Covariate Balance before and after PSM

	Type	Diff. Un	Diff. Adj	M. Threshold
Distance	Distance	0.914 212	0.121 984
GENDER	Binary	0.187 867	0.002 331	Balanced, <0.1
WKWHITE	Binary	0.071 896	0.006 993	Balanced, <0.1
WKSESL	Contin.	−0.306 68	−0.010 52	Balanced, <0.1
RIRT	Contin.	−0.8759	−0.085 21	Balanced, <0.1
MIRT	Contin.	−0.800 92	−0.107 08	Not Balanced, >0.1
S2KPUPRI	Binary	0.181 708	−0.002 33	Balanced, <0.1
P1EXPECT	Contin.	−0.2932	−0.030 87	Balanced, <0.1
P1FIRKDG	Binary	−0.073 43	−0.006 99	Balanced, <0.1
P1AGEENT	Contin.	0.076 272	−0.036 25	Balanced, <0.1
apprchT1	Contin.	−0.643 79	−0.0074	Balanced, <0.1
P1HSEVER	Binary	0.063 855	0	Balanced, <0.1
chg14	Binary	0.004 93	−0.023 31	Balanced, <0.1
avg_RIRT	Contin.	−0.259 18	−0.058	Balanced, <0.1
avg_MIRT	Contin.	−0.197 45	−0.038 75	Balanced, <0.1
avg_SES	Contin.	−0.173 65	−0.011 57	Balanced, <0.1
avg_apprchT1	Contin.	−0.157 44	0.004 81	Balanced, <0.1
S2KMINOR	Contin.	−0.223 03	−0.063 83	Balanced, <0.1
P1FSTAMP	Binary	0.038 837	−0.011 66	Balanced, <0.1
ONEPARENT	Binary	0.050 212	−0.016 32	Balanced, <0.1
STEPPARENT	Binary	0.011 202	−0.016 32	Balanced, <0.1
P1NUMSIB	Contin.	0.146 781	0.002 01	Balanced, <0.1
P1HMAFB	Contin.	−0.262 22	−0.018 13	Balanced, <0.1
WKCAREPK	Binary	−0.026 93	0.051 282	Balanced, <0.1
P1EARLY	Contin.	0.139 938	−0.029 11	Balanced, <0.1
wt_ounces	Contin.	−0.101 39	0.037 886	Balanced, <0.1
C1FMOTOR	Contin.	−0.564 97	−0.050 86	Balanced, <0.1
C1GMOTOR	Contin.	−0.348 13	−0.022 69	Balanced, <0.1
P1HSCALE	Contin.	0.1138 33	0.036 702	Balanced, <0.1
P1SADLON	Contin.	0.083 572	0.020 379	Balanced, <0.1
P1IMPULS	Contin.	0.337 995	−0.017 83	Balanced, <0.1
P1ATTENI	Contin.	0.601 225	−0.012 81	Balanced, <0.1
P1SOLVE	Contin.	0.554 679	0.060 289	Balanced, <0.1
P1PRONOU	Contin.	0.707 568	0.037 271	Balanced, <0.1
P1DISABL	Binary	0.270 212	−0.011 66	Balanced, <0.1

Appendix B: Replication R Code Used in Analysis

B.1. Propensity Score-Based Estimation

library(ggplot2)

library(MatchIt)

library(lmtest)

library(sandwich)

library(cobalt)

data <- eclsk_c

# Define variable categories

demographic <- c("GENDER", "WKWHITE", "WKSESL")

academic <- c("RIRT", "MIRT", "S2KPUPRI", "P1EXPECT", "P1FIRKDG",

"P1AGEENT", "apprchT1", "P1HSEVER", "chg14")

school <- c("avg_RIRT", "avg_MIRT", "avg_SES", "avg_apprchT1", "S2KMINOR")

family <- c("P1FSTAMP", "ONEPARENT", "STEPPARENT", "P1NUMSIB",

"P1HMAFB", "WKCAREPK")

health <- c("P1EARLY", "wt_ounces", "C1FMOTOR", "C1GMOTOR")

parent <- c("P1HSCALE", "P1SADLON", "P1IMPULS", "P1ATTENI",

"P1SOLVE", "P1PRONOU", "P1DISABL")

covariates <- c(demographic, academic, school, family, health, parent)

treatment <- "F5SPECS"

outcome <- "C6R4MSCL"

# Estimate propensity scores using logistic regression

ps_formula <- as.formula(paste(treatment, "~", paste(covariates, collapse = " + ")))

ps_model <- glm(ps_formula, data = data, family = binomial)

data$pscore <- predict(ps_model, type = "response")

# Visualize propensity score distribution

ggplot(data, aes(x = pscore, fill = as.factor(F5SPECS))) +

geom_histogram(position = "identity", bins = 40, alpha = 0.5) +

labs(

title = "Distribution of Estimated Propensity Scores by Treatment Status",

x = "Propensity Score",

y = "Count",

fill = "Treatment"

) +

scale_fill_manual(values = c("#619CFF", "#F8766D"), labels = c("Untreated", "Treated")) +

theme_minimal()

# Perform nearest neighbor matching

psm_model <- matchit(

formula = ps_formula,

data = data,

method = "nearest",

distance = "logit",

ratio = 1

)

summary(psm_model)

# Extract matched dataset

matched_data <- match.data(psm_model)

# Calculate propensity score ranges

min_treated <- min(data$pscore[data[[treatment]] == 1], na.rm = TRUE)

max_treated <- max(data$pscore[data[[treatment]] == 1], na.rm = TRUE)

min_control <- min(data$pscore[data[[treatment]] == 0], na.rm = TRUE)

max_control <- max(data$pscore[data[[treatment]] == 0], na.rm = TRUE)

# Overlap region

overlap_min <- max(min_treated, min_control)

overlap_max <- min(max_treated, max_control)

cat("Treated range:", round(min_treated, 3), "-", round(max_treated, 3), "\n")

cat("Control range:", round(min_control, 3), "-", round(max_control, 3), "\n")

cat("Common support interval:", round(overlap_min, 3), "-", round(overlap_max, 3), "\n")

# Proportion trimmed after matching

trimmed_prop <- 1 - (nrow(matched_data) / nrow(data))

cat("Proportion trimmed:", round(trimmed_prop * 100, 2), "%\n")

# Visualize covariate balance before and after matching (Love plot)

bal <- bal.tab(psm_model)

love.plot(

bal,

threshold = 0.1,

var.order = "unadjusted",

stat = "mean.diffs",

abs = TRUE,

stars = "raw",

title = "Covariate Balance Before and After PSM"

)

bal_tab <- bal.tab(psm_model, un = TRUE, m.threshold = 0.1)

bal_df <- as.data.frame(bal_tab$Balance)

write.csv(bal_df, "balance_statistics_psm.csv", row.names = TRUE)

# Estimate ATE using matched sample (linear regression)

ate_psm <- lm(C6R4MSCL ~ F5SPECS, data = matched_data)

# Compute robust standard errors and 95% confidence interval

robust_se_psm <- coeftest(ate_psm, vcov = vcovHC(ate_psm, type = "HC1"))

ate_coef_psm <- coef(ate_psm)["F5SPECS"]

se_coef_psm <- robust_se_psm["F5SPECS", "Std. Error"]

ci_lower_psm <- ate_coef_psm - 1.96 * se_coef_psm

ci_upper_psm <- ate_coef_psm + 1.96 * se_coef_psm

# Output ATE results

ate_psm_result <- data.frame(

Method = "PSM",

ATE_Estimate = round(ate_coef_psm, 3),

Std_Error = round(se_coef_psm, 3),

CI_Lower = round(ci_lower_psm, 3),

CI_Upper = round(ci_upper_psm, 3)

)

print(ate_psm_result)

# Cluster at the matched pair level (subclass)

cluster_vcov <- vcovCL(ate_psm, cluster = ~subclass, data = matched_data, type = "HC1")

cluster_test <- coeftest(ate_psm, vcov = cluster_vcov)

# Extract cluster-robust ATE results

ate_cluster_coef <- cluster_test["F5SPECS", "Estimate"]

ate_cluster_se <- cluster_test["F5SPECS", "Std. Error"]

ci_lower_cluster <- ate_cluster_coef - 1.96 * ate_cluster_se

ci_upper_cluster <- ate_cluster_coef + 1.96 * ate_cluster_se

ate_psm_cluster_result <- data.frame(

Method = "PSM (Cluster-Robust SE)",

ATE_Estimate = round(ate_cluster_coef, 3),

Std_Error = round(ate_cluster_se, 3),

CI_Lower = round(ci_lower_cluster, 3),

CI_Upper = round(ci_upper_cluster, 3)

)

print(ate_psm_cluster_result)

B.2. Machine Learning-Based Estimation

B.2.1. Bayesian Additive Regression Trees (BART.R)

library(BART)

library(ggplot2)

# Define covariates, treatment, outcome

X_bart <- data[, covariates] # Covariate matrix

W_bart <- data[[treatment]] # Treatment indicator

Y_bart <- data[[outcome]] # Outcome variable

# Split data by treatment status

X_treated <- X_bart[W_bart == 1, ]

Y_treated <- Y_bart[W_bart == 1]

X_control <- X_bart[W_bart == 0, ]

Y_control <- Y_bart[W_bart == 0]

# Fit BART models separately for treated and control groups

set.seed(2025)

bart_model_treated <- gbart(x.train = X_treated, y.train = Y_treated, x.test = X_bart)

set.seed(2025)

bart_model_control <- gbart(x.train = X_control, y.train = Y_control, x.test = X_bart)

# Extract predicted outcomes

mu1_hat <- colMeans(bart_model_treated$yhat.test) # Predicted outcome if treated

mu0_hat <- colMeans(bart_model_control$yhat.test) # Predicted outcome if untreated

# Estimate ATE and 95% CI

ite_bart_values <- mu1_hat - mu0_hat

ate_bart_estimate <- mean(ite_bart_values)

ate_bart_se <- sd(ite_bart_values) / sqrt(length(ite_bart_values))

ci_lower_bart <- ate_bart_estimate - 1.96 * ate_bart_se

ci_upper_bart <- ate_bart_estimate + 1.96 * ate_bart_se

ate_bart_result <- data.frame(

Method = "BART",

ATE_Estimate = round(ate_bart_estimate, 3),

Std_Error = round(ate_bart_se, 3),

CI_Lower = round(ci_lower_bart, 3),

CI_Upper = round(ci_upper_bart, 3)

)

print(ate_bart_result)

# Plot distribution of ITEs

ite_bart_df <- data.frame(ITE = ite_bart_values)

ggplot(ite_bart_df, aes(x = ITE)) +

geom_histogram(bins = 40, fill = "firebrick", alpha = 0.7, color = "black") +

labs(title = "Distribution of Individual Treatment Effects (ITE) from BART",

x = "Individual Treatment Effect",

y = "Frequency") +

theme_minimal()

ite_bart_summary <- summary(ite_bart_values)

ite_bart_sd <- sd(ite_bart_values)

print(ite_bart_summary)

print(ite_bart_sd)

B.2.2. Causal Forest (CF.R)

library(grf)

library(ggplot2)

set.seed(2025)

# Define inputs for Causal Forest

X_cf <- data[, covariates] # Covariate matrix

W_cf <- data[[treatment]] # Treatment indicator

Y_cf <- data[[outcome]] # Outcome variable

# Fit Causal Forest model

cf_model <- causal_forest(X = X_cf, Y = Y_cf, W = W_cf)

# Estimate ATE and 95% confidence interval

ate_cf_est <- average_treatment_effect(cf_model, target.sample = "all")

ate_cf_result <- data.frame(

Method = "Causal Forest",

ATE_Estimate = round(ate_cf_est[["estimate"]], 3),

Std_Error = round(ate_cf_est[["std.err"]], 3),

CI_Lower = round(ate_cf_est[["estimate"]] - 1.96 * ate_cf_est[["std.err"]], 3),

CI_Upper = round(ate_cf_est[["estimate"]] + 1.96 * ate_cf_est[["std.err"]], 3)

)

print(ate_cf_result)

# Predict Individual Treatment Effects (ITE) from the fitted Causal Forest model

cf_ite_pred <- predict(cf_model, estimate.variance = TRUE)

# Plot the distribution of ITE

ite_cf_df <- data.frame(ITE = ite_cf_values)

ggplot(ite_cf_df, aes(x = ITE)) +

geom_histogram(bins = 40, fill = "steelblue", alpha = 0.7, color = "black") +

labs(title = "Distribution of Individual Treatment Effects (ITE) from Causal Forest",

x = "Individual Treatment Effect",

y = "Frequency") +

theme_minimal()

# Extract ITE predictions and variances

ite_cf_values <- cf_ite_pred$predictions

ite_cf_variances <- cf_ite_pred$variance.estimates

ite_summary <- summary(ite_cf_values)

ite_sd <- sd(ite_cf_values)

print(ite_summary)

print(ite_sd)

B.2.3. Double Machine Learning (DML.R)

library(DoubleML)

library(mlr3)

library(mlr3learners)

library(data.table)

set.seed(2025)

# Prepare data

dml_data <- data[, c(treatment, outcome, covariates)]

setDT(dml_data)

dml_data_obj <- DoubleMLData$new(data = dml_data,

y_col = outcome,

d_cols = treatment)

# Define learners

learner_l <- lrn("regr.ranger", num.trees = 500, max.depth = 5) # For outcome

learner_m <- lrn("classif.ranger", predict_type = "prob", num.trees = 500, max.depth = 5) # For treatment

# Instantiate DML model

dml_plr <- DoubleMLPLR$new(data = dml_data_obj,

ml_l = learner_l,

ml_m = learner_m,

n_folds = 5,

score = "partialling out")

dml_plr$fit()

ate_dml_result <- data.frame(

Method = "DML",

ATE_Estimate = round(dml_plr$coef, 3),

Std_Error = round(dml_plr$se, 3),

CI_Lower = round(dml_plr$coef - 1.96 * dml_plr$se, 3),

CI_Upper = round(dml_plr$coef + 1.96 * dml_plr$se, 3)

)

print(ate_dml_result)

# Prepare two new datasets: one with all treatment = 1, one with all treatment = 0

dml_data_treat1 <- copy(dml_data)

dml_data_treat1[[treatment]] <- 1

dml_data_treat0 <- copy(dml_data)

dml_data_treat0[[treatment]] <- 0

# Train outcome model learner on original dataset

task_dml <- TaskRegr$new(id = "outcome_dml", backend = dml_data, target = outcome)

learner_l$train(task_dml)

# Predict potential outcomes under treatment and control

y_hat_dml_1 <- learner_l$predict_newdata(newdata = dml_data_treat1)$response

y_hat_dml_0 <- learner_l$predict_newdata(newdata = dml_data_treat0)$response

# Calculate ITEs and plot distribution

ite_dml_values <- y_hat_dml_1 - y_hat_dml_0

ite_dml_df <- data.frame(ITE = ite_dml_values)

ggplot(ite_dml_df, aes(x = ITE)) +

geom_histogram(bins = 40, fill = "darkorange", alpha = 0.7, color = "black") +

labs(title = "Distribution of Individual Treatment Effects (ITE) from DML",

x = "Individual Treatment Effect",

y = "Frequency") +

theme_minimal()

ite_dml_summary <- summary(ite_dml_values)

ite_dml_sd <- sd(ite_dml_values)

print(ite_dml_summary)

print(ite_dml_sd)

B.3. Additional Analyses

library(randomForest)

set.seed(2025) # Ensure reproducibility

# Fit a random forest model to predict fifth-grade math scores

rf_model <- randomForest(

C6R4MSCL ~ .,

data = data[, c(covariates, "C6R4MSCL")],

importance = TRUE

)

# Extract variable importance (type = 1 corresponds to %IncMSE)

importance_df <- as.data.frame(importance(rf_model, type = 1))

importance_df$Variable <- rownames(importance_df)

# Sort by %IncMSE in descending order

importance_df <- importance_df[order(-importance_df[["%IncMSE"]]), ]

# Select top 10 most important variables

top_vars <- head(importance_df, 10)

# Plot variable importance

ggplot(top_vars, aes(x = reorder(Variable, `%IncMSE`), y = `%IncMSE`, fill = `%IncMSE`)) +

geom_col(width = 0.7, show.legend = FALSE) + # Bars

geom_text(aes(label = round(`%IncMSE`, 1)), hjust = -0.1, size = 4.2) + # Value labels

coord_flip() + # Horizontal bars

scale_fill_gradient(low = "#dceaf5", high = "#1f78b4") + # Soft blue gradient

labs(

title = "Top 10 Most Important Variables for Predicting Math Scores",

x = NULL,

y = "Increase in MSE (%)"

) +

theme_minimal(base_size = 10) +

theme(

plot.title = element_text(hjust = 0.5, face = "bold"),

axis.text.y = element_text(face = "bold"),

panel.grid.major.y = element_blank()

) +

scale_y_continuous(expand = expansion(mult = c(0, 0.1)))

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Rubin, D.B. (1974) Estimating Causal Effects of Treatments in Randomized and Nonrandomized Studies. Journal of Educational Psychology, 66, 688-701.[CrossRef]
[2]	Holland, P.W. (1986) Statistics and Causal Inference. Journal of the American Statistical Association, 81, 945-960.[CrossRef]
[3]	Rubin, D.B. (1986) Comment: Which Ifs Have Causal Answers. Journal of the American Statistical Association, 81, 961-962.[CrossRef]
[4]	Imbens, G.W. and Rubin, D.B. (2015) Causal Inference for Statistics, Social, and Biomedical Sciences. Cambridge University Press.[CrossRef]
[5]	Spirtes, P., Glymour, C. and Scheines, R. (2001) Causation, Prediction, and Search. 2nd Edition, The MIT Press.[CrossRef]
[6]	Belloni, A., Chernozhukov, V., Fernandez-Val, I. and Hansen, C. (2017) Program Evaluation and Causal Inference with High-Dimensional Data. Econometrica, 85, 233-298.[CrossRef]
[7]	Fan, J., Han, F. and Liu, H. (2014) Challenges of Big Data Analysis. National Science Review, 1, 293-314.[CrossRef] [PubMed]
[8]	Filzmoser, P. and Nordhausen, K. (2020) Robust Linear Regression for High-Dimensional Data: An Overview. WIREs Computational Statistics, 13, e1524.[CrossRef]
[9]	Kleinberg, J., Ludwig, J., Mullainathan, S. and Obermeyer, Z. (2015) Prediction Policy Problems. American Economic Review, 105, 491-495.[CrossRef] [PubMed]
[10]	Rose, S. (2016) A Machine Learning Framework for Plan Payment Risk Adjustment. Health Services Research, 51, 2358-2374.[CrossRef] [PubMed]
[11]	Ryll, L. and Seidens, S. (2019) Evaluating the Performance of Machine Learning Algorithms in Financial Market Forecasting: A Comprehensive Survey. arXiv: 1906.07786.
[12]	McConnell, K.J. and Lindner, S. (2019) Estimating Treatment Effects with Machine Learning. Health Services Research, 54, 1273-1282.[CrossRef] [PubMed]
[13]	Opper, I.M. (2020) Using Machine Learning Techniques to Improving Average Treatment Effect Estimates in Small-Scale Randomized Controlled Trials. RAND.
[14]	Lake, B.M., Ullman, T.D., Tenenbaum, J.B. and Gershman, S.J. (2016) Building Machines That Learn and Think Like People. Behavioral and Brain Sciences, 40, e253.[CrossRef] [PubMed]
[15]	Sagawa, S., Raghunathan, A., Koh, P.W. and Liang, P. (2020) An Investigation of Why Overparameterization Exacerbates Spurious Correlations. International Conference on Machine Learning, 13-18 July 2020, 8346-8356.
[16]	Salahuddin, Z., Woodruff, H.C., Chatterjee, A. and Lambin, P. (2022) Transparency of Deep Neural Networks for Medical Image Analysis: A Review of Interpretability Methods. Computers in Biology and Medicine, 140, Article ID: 105111.[CrossRef] [PubMed]
[17]	Pearl, J. (2009) Causality. 2nd Edition, Cambridge University Press.[CrossRef]
[18]	Feuerriegel, S., Frauen, D., Melnychuk, V., Schweisthal, J., Hess, K., Curth, A., et al. (2024) Causal Machine Learning for Predicting Treatment Outcomes. Nature Medicine, 30, 958-968.[CrossRef] [PubMed]
[19]	Hill, J.L. (2011) Bayesian Nonparametric Modeling for Causal Inference. Journal of Computational and Graphical Statistics, 20, 217-240.[CrossRef]
[20]	Stanovich, K.E. (2009) Matthew Effects in Reading: Some Consequences of Individual Differences in the Acquisition of Literacy. Journal of Education, 189, 23-55.[CrossRef]
[21]	Sanchez, P., Voisey, J.P., Xia, T., Watson, H.I., O’Neil, A.Q. and Tsaftaris, S.A. (2022) Causal Machine Learning for Healthcare and Precision Medicine. Royal Society Open Science, 9, Article ID: 220638.[CrossRef] [PubMed]
[22]	Kuzmanovic, M., Frauen, D., Hatt, T. and Feuerriegel, S. (2024) Causal Machine Learning for Cost-Effective Allocation of Development Aid. Proceedings of the 30th ACM SIGKDD Conference on Knowledge Discovery and Data Mining, Barcelona, 25-29 August 2024, 5283-5294.[CrossRef]
[23]	Hwang, N. (2024) Label to Help: The Effects of Special Education on Student Outcomes. EdWorkingPaper No. 24-939. Annenberg Institute for School Reform at Brown University.
[24]	Hornby, G. and Kauffman, J.M. (2024) Inclusive Education, Intellectual Disabilities and the Demise of Full Inclusion. Journal of Intelligence, 12, Article 20.[CrossRef] [PubMed]
[25]	Stark, K., Koslouski, J.B., Vadhan, J. and Vega, M. (2024) The Future Is Inclusive: An Invitation for Interdisciplinary Collaboration between Social Emotional Learning and Special Education Researchers. Social and Emotional Learning: Research, Practice, and Policy, 3, Article ID: 100043.[CrossRef]
[26]	Dell’Anna, S., Entrich, S.R. and Banks, J. (2025) Assessing the Outcomes of Students with Special Educational Needs in Inclusive Education: A Comparative Study of Germany, Ireland, and Italy. PROSPECTS.[CrossRef]
[27]	Rosenbaum, P.R. and Rubin, D.B. (1983) The Central Role of the Propensity Score in Observational Studies for Causal Effects. Biometrika, 70, 41-55.[CrossRef]
[28]	Chipman, H.A., George, E.I. and McCulloch, R.E. (2010) BART: Bayesian Additive Regression Trees. The Annals of Applied Statistics, 4, 266-298.[CrossRef]
[29]	Wager, S. and Athey, S. (2018) Estimation and Inference of Heterogeneous Treatment Effects Using Random Forests. Journal of the American Statistical Association, 113, 1228-1242.[CrossRef]
[30]	Chernozhukov, V., Chetverikov, D., Demirer, M., Duflo, E., Hansen, C., Newey, W. and Robins, J. (2018) Double/Debiased Machine Learning for Treatment and Structural Parameters. The Econometrics Journal, 21, C1-C68.
[31]	Winship, C. and Radbill, L. (1994) Sampling Weights and Regression Analysis. Sociological Methods & Research, 23, 230-257.[CrossRef]
[32]	Solon, G., Haider, S.J. and Wooldridge, J.M. (2015) What Are We Weighting For? Journal of Human Resources, 50, 301-316.[CrossRef]
[33]	Breiman, L. (2001) Random Forests. Machine Learning, 45, 5-32.[CrossRef]