When Prediction Error Is Not Enough: Evaluating Nuisance-Function Prediction for Causal Estimation

Abstract

Prediction error is widely used to evaluate nuisance-function estimators in causal inference, but its relationship with causal estimator performance may differ across performance measures. We studied this question in a partially linear model using Monte Carlo simulations. We compared ordinary least squares (OLS), generalized additive models (GAMs), XGBoost, and Double Machine Learning with XGBoost (DML-XGBoost), evaluating nuisance-function prediction error, bias, RMSE, and 95% confidence interval coverage. We also examined a simple joint-error measure based on the absolute cross-product of estimation errors from the exposure and outcome nuisance functions. Across the simulated settings, XGBoost had the lowest RMSE among the non-oracle methods, while DML-XGBoost generally provided better confidence interval coverage. Prediction error did not consistently track causal bias across methods and settings, and the method with the best point-estimation performance did not necessarily have the best confidence interval coverage. The joint-error measure was only weakly associated with causal bias and did not provide a useful standalone measure of causal performance. These results suggest that prediction error is useful for assessing nuisance-function estimation, but it should not be treated as a direct measure of the quality of the resulting causal estimator.

Keywords: causal inference; nuisance functions; double machine learning; causal estimation

1 Introduction

Causal effect estimation with observational data often requires estimating nuisance functions that describe the relationships between covariates, the exposure, and the outcome. Flexible regression and machine-learning methods allow these functions to be estimated without imposing restrictive parametric assumptions, and they are now widely used in causal inference [8, 9]. Double Machine Learning (DML), for example, combines flexible nuisance-function estimation with Neyman-orthogonal scores and cross-fitting to reduce the impact of nuisance estimation error on the target causal parameter [3, 5, 12, 4, 1]. In practice, however, nuisance models are often evaluated using the same prediction measures used in ordinary supervised learning. This raises a simple question: does better prediction of the nuisance functions also mean better causal estimation?

Prediction error is a natural way to evaluate nuisance models because estimating a nuisance function is itself a prediction problem. Cross-validation error and mean squared error are commonly used to compare learners and tune their hyperparameters [7]. Recent studies have also considered whether prediction performance can help with causal analysis. Bach et al. [2] studied learner choice, hyperparameter tuning, and data-splitting decisions in Double Machine Learning. They found that predictive performance can be useful for model selection, but its relationship with causal performance depends on the learner and the causal setting. Doutreligne and Varoquaux [6] compared standard predictive criteria, including mean squared error, with criteria based more directly on causal objectives using simulations and several health-care datasets. These studies focus mainly on using prediction performance to choose models or learners. We instead ask whether the prediction error of a fitted nuisance function can be used to anticipate how well the resulting causal estimator performs. This question is important because good predictive performance does not necessarily imply good performance for a causal estimand [11, 10].

The issue is more complicated when a causal estimator depends on more than one nuisance function. In the partially linear model with a continuous exposure, both the exposure regression and the outcome regression enter the DML score [3, 4]. Their errors can therefore affect the causal estimator together. Orthogonality reduces the leading first-order effect of nuisance estimation error, but it does not remove finite-sample differences caused by nuisance estimation [3, 12, 1]. We therefore examine both the prediction error of each nuisance function and a simple measure of their joint error. We then ask whether these measures are related to finite-sample causal performance. The prediction error is measured for the nuisance functions, while bias, RMSE, and confidence interval coverage are measured for the causal estimator. Thus, the two sets of measures capture different parts of the analysis.

We study these questions using Monte Carlo simulation in a partially linear causal model with continuous exposure. We compare ordinary least squares (OLS), generalized additive models (GAMs), XGBoost, and Double Machine Learning with XGBoost (DML-XGBoost). We evaluate nuisance-function prediction error, causal bias, RMSE, and 95% confidence interval coverage. We also examine whether a simple measure of joint nuisance-function error is related to finite-sample causal bias. Finally, we consider a clustered-data setting with cluster-level cross-fitting to determine whether the main findings change when observations within clusters are dependent.

2 Methods

2.1 Data-Generating Framework

We conducted a Monte Carlo simulation study to examine whether nuisance-function prediction accuracy was associated with causal effect estimation performance. The simulation compared prediction error, causal estimation accuracy, and inferential performance across several nuisance- function estimation methods and two dependence settings.

Let $X_{i}$ denote a continuous exposure, $Y_{i}$ the outcome, and $C_{i}=(C_{i1},\ldots,C_{ip})^{\top}$ a vector of observed baseline covariates. We generated data from the partially linear model

$$X_{i}=m_{0}(C_{i})+\epsilon_{Xi},$$

$$Y_{i}=\tau X_{i}+g_{0}(C_{i})+\epsilon_{Yi},$$

where

$$m_{0}(C_{i})=E(X_{i}\mid C_{i})$$

is the exposure nuisance function and $g_{0}(C_{i})$ is the outcome baseline regression function. The causal effect $\tau$ was constant across individuals. We assumed

$$E(\epsilon_{Xi}\mid C_{i})=0,\qquad E(\epsilon_{Yi}\mid X_{i},C_{i})=0.$$

For the outcome regression entering the partially linear orthogonal score, the corresponding nuisance function is

$$h_{0}(C_{i})=E(Y_{i}\mid C_{i})=\tau m_{0}(C_{i})+g_{0}(C_{i}).$$

Thus, the two nuisance functions used in the causal estimator were $m_{0}(C)$ and $h_{0}(C)$. The true functions were known by construction, allowing prediction error to be evaluated directly rather than through noisy observed outcomes.

2.2 Nuisance-Function Error and Causal Estimation

For fitted nuisance functions $\widehat{m}(C)$ and $\widehat{h}(C)$, where $m_{0}(C)=E(X\mid C)$ and $h_{0}(C)=E(Y\mid C)=\tau m_{0}(C)+g_{0}(C)$, we defined the prediction errors

$$\delta_{m}(C)=\widehat{m}(C)-m_{0}(C),\qquad\delta_{h}(C)=\widehat{h}(C)-h_{0}(C).$$

We summarized nuisance-function prediction error using the mean squared prediction errors

$$\mathcal{E}_{m}=E\{\delta_{m}(C)^{2}\},\qquad\mathcal{E}_{h}=E\{\delta_{h}(C)^{2}\},$$

The evaluation sample was not used for fitting any nuisance model or for estimating the causal effect. Thus, the reported nuisance-function errors represent out-of-sample discrepancies from the known data-generating functions rather than training-set prediction errors. These quantities measure the average magnitude of prediction error over the covariate distribution but contain no information about where that error occurs.

We additionally considered a simple descriptive measure of joint nuisance-function error,

$$\mathcal{D}_{\mathrm{joint}}=\left|E\{\delta_{m}(C)\delta_{h}(C)\}\right|.$$

The signed quantity $E\{\delta_{m}(C)\delta_{h}(C)\}$ describes the direction of the average alignment between the two nuisance-function errors, whereas $\mathcal{D}_{\mathrm{joint}}$ discards this direction and summarizes its absolute magnitude. A value near zero may therefore arise either because the two errors are small or because positive and negative contributions cancel across the covariate distribution.

We compared three nuisance-function learners—ordinary least squares (OLS), generalized additive models (GAMs), and XGBoost—using the same residualized treatment-effect estimator. OLS provided a parametric benchmark, GAMs allowed nonlinear additive relationships, and XGBoost allowed nonlinearities and interactions. For each learner, separate models were fitted for the exposure nuisance function $m(C)$, and the outcome nuisance function $h(C)$, and the resulting residualized exposure and outcome were used to estimate $\tau$. We additionally considered Double Machine Learning with XGBoost (DML-XGBoost), in which the same XGBoost learner was used within a cross-fitted orthogonal estimation procedure. Thus, the comparison between XGBoost and DML-XGBoost reflects not only nuisance-function learning but also the use of cross-fitting and orthogonalization.

For OLS, GAM, and XGBoost, nuisance functions were estimated using the full training sample. For DML-XGBoost, nuisance functions were estimated using the training observations outside the corresponding cross-fitting fold. For all methods, prediction error was evaluated on the same independent evaluation sample. For DML-XGBoost, prediction error was calculated by applying each cross-fitting-specific nuisance estimate to the independent evaluation sample and averaging the resulting squared prediction errors across cross-fitting fits. For observation $i$ in fold $k$, the DML nuisance predictions

$$\widehat{m}^{(-k)}(C_{i})\quad\text{and}\quad\widehat{h}^{(-k)}(C_{i})$$

were obtained using observations outside fold $k$. The causal effect was then estimated from the cross-fitted residualized variables

$$\widetilde{X}_{i}=X_{i}-\widehat{m}^{(-k(i))}(C_{i}),\qquad\widetilde{Y}_{i}=Y_{i}-\widehat{h}^{(-k(i))}(C_{i}),$$

using

$$\widehat{\tau}=\frac{\sum_{i=1}^{n}\widetilde{X}_{i}\widetilde{Y}_{i}}{\sum_{i=1}^{n}\widetilde{X}_{i}^{2}}.$$

Cross-fitting ensured that each observation was evaluated using nuisance estimates obtained without using that observation for nuisance-model fitting. In clustered settings, all observations from the same cluster were assigned to the same fold.

We additionally considered a clustered data-generating setting to evaluate the sensitivity of nuisance-function estimation and causal inference under within-cluster dependence. Let $j(i)$ denote the cluster containing individual $i$. We generated

$$X_{i}=m_{0}(C_{i})+u^{X}_{j(i)}+\epsilon_{Xi},$$

and

$$Y_{i}=\tau X_{i}+g_{0}(C_{i})+u^{Y}_{j(i)}+\epsilon_{Yi},$$

where $u^{X}_{j}$ and $u^{Y}_{j}$ are cluster-level random effects shared by individuals within cluster $j$. These terms induce within-cluster dependence in both the exposure and outcome while preserving a constant causal effect $\tau$. The cluster-level random effects were generated independently of the covariates and individual-level errors and had mean zero.

For the clustered scenario, the training sample contained 1,000 observations organized into clusters of size 10. Cluster-level random effects had standard deviation 0.8 for both the exposure and outcome models, and the individual-level error standard deviations were 1.0. Five-fold cross-fitting was performed at the cluster level, such that all observations from a given cluster were assigned to the same fold. Causal inference used cluster-level score contributions for variance estimation. The clustered scenario was treated as a secondary dependence-structure analysis rather than a systematic evaluation of clustering strength.

2.3 Simulation Design and Performance Evaluation

We considered two simulation scenarios. Scenario 1 represented the baseline setting with independent observations. Scenario 2 represented the clustered setting, in which observations were nested within clusters and cluster-level random effects induced within-cluster dependence. For each scenario, the covariate vector $C_{i}$ was generated from the same distribution used to define the underlying nuisance functions, and the causal effect was set to a constant value $\tau$. The nonlinear components of $m_{0}(C)$ and $g_{0}(C)$ were fixed across Monte Carlo replicates so that differences across methods reflected estimation performance rather than changes in the data-generating mechanism. The simulation was organized around two primary questions. First, we examined whether nuisance-function prediction error was informative about causal estimation performance. Second, we examined whether a simple summary of joint nuisance-function error provided information about finite-sample causal bias beyond the marginal prediction errors.

We additionally considered variation in fitted nuisance functions across repeated training samples as a complementary property of nuisance-function estimation, but did not evaluate it as a separate simulation endpoint. The clustered scenario additionally evaluated these relationships under within-cluster dependence. The training sample size was $n=1000$ in both scenarios. Within each replicate, the training data were used for nuisance-function estimation and causal effect estimation, whereas an independent evaluation sample was used to assess nuisance-function prediction error against the known true nuisance functions. The same simulated dataset was used across competing methods within each replicate to facilitate comparison of finite-sample performance. For each method and simulation setting, we recorded prediction error, the joint nuisance-error measure, bias, RMSE, empirical standard error, confidence interval coverage, and confidence interval width.Sample-to-sample variation in fitted nuisance functions was considered as a complementary conceptual property of nuisance-function estimation, but was not evaluated as a separate simulation endpoint. For a causal estimator $\widehat{\tau}$, bias and RMSE were calculated as

$$\mathrm{Bias}=E(\widehat{\tau})-\tau$$

and

$$\mathrm{RMSE}=\left[E\{(\widehat{\tau}-\tau)^{2}\}\right]^{1/2}.$$

Confidence interval coverage was defined as the proportion of Monte Carlo replicates in which the confidence interval contained the true value of $\tau$.

2.4 Inference

For the partially linear estimator, inference was based on the estimated influence function associated with the orthogonal moment. Let

$$\widetilde{X}_{i}=X_{i}-\widehat{m}_{i},\qquad\widetilde{Y}_{i}=Y_{i}-\widehat{h}_{i},$$

and define

$$\widehat{\psi}_{i}=\widetilde{X}_{i}\left\{\widetilde{Y}_{i}-\widehat{\tau}\widetilde{X}_{i}\right\}.$$

The empirical variance of the influence function was combined with the derivative of the estimating equation to obtain the standard error of $\widehat{\tau}$. Wald-type 95% confidence intervals were constructed as

$$\widehat{\tau}\pm 1.96\,\widehat{\mathrm{SE}}(\widehat{\tau}).$$

For each simulation setting, we reported Monte Carlo bias, RMSE, empirical standard error, estimated standard error, and confidence interval coverage. The empirical standard error was calculated as the standard deviation of the estimated causal effects across Monte Carlo replicates. Coverage was calculated as the proportion of replicates in which the nominal 95% confidence interval contained the true causal effect. In the clustered setting, the influence-function variance was estimated using cluster-level score contributions, and the same cluster structure was used for both cross-fitting and variance estimation. The primary analysis compared marginal nuisance-function prediction error with causal estimation performance using bias, RMSE, and confidence interval coverage. We also compared method rankings based on prediction error with rankings based on causal performance.

3 Results

Table 1: Causal estimation performance across simulation scenarios.

Scenario	Method	Bias	RMSE	Empirical SE	Model SE	Coverage
1	Oracle	0.0004	0.0101	0.0101	0.0101	0.950
1	OLS	0.1092	0.1099	0.1099	0.0157	0.000
1	GAM	0.0630	0.0642	0.0642	0.0150	0.000
1	XGBoost	-0.0050	0.0143	0.0144	0.0130	0.856
1	DML-XGBoost	-0.0057	0.0372	0.0370	0.0340	0.928
2	Oracle	-0.0005	0.0181	0.0181	0.0181	0.954
2	OLS	0.0680	0.0701	0.0701	0.0157	0.018
2	GAM	0.0382	0.0417	0.0417	0.0150	0.300
2	XGBoost	-0.0104	0.0208	0.0208	0.0180	0.838
2	DML-XGBoost	-0.0105	0.0522	0.0511	0.0450	0.908

Note. Bias is the Monte Carlo mean of $\widehat{\tau}-\tau$; RMSE is the root mean squared error across replicates. Empirical SE is the standard deviation of $\widehat{\tau}$ across replicates, and Model SE is the mean estimated standard error. Coverage is the proportion of nominal 95% confidence intervals containing the true causal effect. Scenario 2 incorporates cluster-level dependence with cluster-level cross-fitting and cluster-robust variance estimation.

Across the two simulation settings, XGBoost had the lowest RMSE among the non-oracle methods, whereas DML-XGBoost provided better confidence interval coverage. In Scenario 1, XGBoost had an RMSE of 0.0143 and bias of -0.0050, compared with an RMSE of 0.0372 and bias of -0.0057 for DML-XGBoost. In Scenario 2, the corresponding RMSE values were 0.0208 and 0.0522, with biases of 0.0104 and 0.0105, respectively. XGBoost had confidence interval coverage of 0.856 and 0.838 in Scenarios 1 and 2, respectively, whereas DML-XGBoost had coverage of 0.928 and 0.908. The oracle estimator had coverage of 0.950 and 0.954, with biases of 0.0004 and -0.0005, respectively. Thus, the method with the best point-estimation performance was not the method with the best inferential calibration.

Prediction error did not give the same ranking as causal performance. In Scenario 1, XGBoost had lower prediction error for both nuisance functions than DML-XGBoost, with mean $\mathcal{E}_{m}$ of 0.169 and mean $\mathcal{E}_{h}$ of 0.412, compared with 0.190 and 0.467 for DML-XGBoost. In Scenario 2, the corresponding mean $\mathcal{E}_{m}$. Thus, lower nuisance-function prediction error was not sufficient to identify the method with better inferential performance. Thus, nuisance-function prediction error correctly identified the stronger point-estimation method in these simulations, but it did not identify the method with the best inferential calibration.

The joint nuisance-error measure $\mathcal{D}_{\mathrm{joint}}=|E\{\delta_{m}(C)\delta_{h}(C)\}|$ showed weak associations with absolute causal bias. Among DML-XGBoost replicates, the correlation between $\mathcal{D}_{\mathrm{joint}}$ and absolute bias was -0.111 in Scenario 1 and -0.065 in Scenario 2. Across the practical learners, the correlations in Scenario 2 ranged from approximately -0.10 to 0.11 (Figure 2). These results do not suggest a strong relationship between this simple joint-error measure and finite-sample causal bias. We therefore view $\mathcal{D}_{\mathrm{joint}}$ as a descriptive diagnostic rather than as a measure for ranking or selecting causal estimators. Rather than serving as a standalone predictor or estimator-selection criterion, the measure is better interpreted as a descriptive diagnostic of the alignment of errors in the two nuisance functions. In the clustered setting, the main pattern was similar. The oracle estimator remained approximately unbiased with close-to-nominal coverage, while the practical learners differed in both point-estimation accuracy and confidence interval coverage. Because only one cluster size and one set of variance components were considered, this analysis was intended as a sensitivity analysis rather than a systematic study of clustering strength.

Figure 1: Relationship between nuisance-function prediction error and absolute causal bias across Monte Carlo replicates. The strength and direction of the associations varied across learners and simulation settings, indicating that prediction error alone did not fully characterize finite-sample causal performance.

Figure 2: Relationship between the joint nuisance-error measure $\mathcal{D}_{\mathrm{joint}}$ and absolute causal bias across Monte Carlo replicates. The weak associations suggest that this simple joint-error measure does not provide a reliable standalone summary of finite-sample causal bias.

4 Discussion

This simulation study shows that nuisance-function prediction performance and causal inferential performance do not always coincide. Across the simulated settings, XGBoost generally achieved the lowest RMSE among the non-oracle methods, whereas DML-XGBoost provided better confidence interval coverage. Thus, point-estimation accuracy and inferential calibration were not the same across methods. These findings suggest that conventional prediction error provides useful information about nuisance-function quality, but should not be used as the sole measure of causal estimation performance. In our simulations, lower prediction error was associated with better point-estimation performance for XGBoost, but it did not correspond to better confidence interval coverage. Prediction error is therefore useful for assessing nuisance-model quality, but it should not be viewed as a standalone surrogate for overall causal estimation quality. More broadly, nuisance-function assessment may involve multiple distinct properties. Prediction error measures accuracy relative to the underlying nuisance function, whereas function-level discrepancies describe variation among fitted functions across repeated training samples. These properties need not rank learning procedures in the same way. A learner could, for example, have low prediction error but substantial sample-to-sample variation, or produce similar fitted functions while retaining systematic prediction error. Thus, prediction accuracy and sample-to-sample variation describe complementary aspects of nuisance-function estimation.

The joint nuisance-error measure $\mathcal{D}_{\mathrm{joint}}$ was only weakly related to absolute causal bias, and the direction of the association varied across methods. In these simulations, the measure was more useful as a descriptive summary of the relationship between the two nuisance-function errors than as a criterion for comparing causal estimators. It should not be interpreted as a complete measure of orthogonal-score error. The clustered analysis gave similar qualitative findings, but its scope was limited. Only one cluster size and one set of variance components were considered, so these results do not establish performance across a wider range of dependence structures. The results should not be interpreted as showing that prediction error is irrelevant for causal inference. Rather, prediction error remains a useful measure of nuisance-function accuracy, particularly for assessing whether a learner is estimating the underlying regression functions adequately. The limitation is that it does not, by itself, determine how nuisance estimation translates into the sampling distribution of the final causal estimator.

Several limitations should be noted. First, the simulation considered a constant treatment effect and a relatively simple partially linear causal structure. Second, the clustered analysis was limited to one cluster size and one set of variance components and therefore does not characterize performance across a broader range of dependence structures. Third, the comparison between DML-XGBoost and the other approaches does not isolate the effect of the nuisance learner because DML-XGBoost combines XGBoost with cross-fitting and an orthogonal estimating procedure. Therefore, differences in causal performance reflect both the nuisance learner and the causal estimation procedure. Finally, the joint-error analysis considered a simple descriptive summary rather than the full structure of the orthogonal-score remainder. Because the simulations considered a constant treatment effect, the relationship between nuisance-function prediction error and causal performance may differ in settings with treatment-effect heterogeneity. Future work could investigate whether variation in fitted nuisance functions across repeated training samples provides information about causal estimator performance beyond conventional prediction error. In particular, this variation could be evaluated across different regions of the covariate distribution and related to downstream bias, RMSE, and confidence interval coverage. Such analyses could also examine whether incorporating sample-to-sample variation into nuisance-function assessment improves the selection of learners for causal estimation.

Data and Code Availability

This study used simulated data generated from the data-generating mechanisms described in the Methods section. No human participant data were used. The simulation and analysis code are publicly available on GitHub at
https://github.com/congca2/nuisance-function-prediction-causal-estimation.

Funding

The author received no specific funding for this work.

Author Contributions

Cong Cao: Conceptualization, Methodology, Formal analysis, Software, Visualization, Writing – original draft, Writing – review and editing.

Competing Interests

The authors declare no competing interests.