Paper deep dive
Analyzing Error Sources in Global Feature Effect Estimation
Timo Heiß, Coco Bögel, Bernd Bischl, Giuseppe Casalicchio
Intelligence
Status: succeeded | Model: google/gemini-3.1-flash-lite-preview | Prompt: intel-v1 | Confidence: 95%
Last extracted: 3/22/2026, 5:17:39 AM
Summary
This paper provides a systematic, estimator-level analysis of error sources in global feature effect estimation, specifically for Partial Dependence (PD) and Accumulated Local Effects (ALE) plots. The authors derive a mean-squared-error (MSE) decomposition that separates model bias, estimation bias, model variance, and estimation variance. They find that while holdout data is theoretically cleaner, training data biases are empirically negligible due to larger sample sizes. They also highlight that ALE is particularly sensitive to sample size and that cross-validation can effectively reduce model variance.
Entities (4)
Relation Signals (3)
Partial Dependence → haserrorcomponents → Mean Squared Error Decomposition
confidence 95% · we derive a mean-squared-error decomposition that separates model bias, estimation bias, model variance, and estimation variance for PD and ALE.
Accumulated Local Effects → haserrorcomponents → Mean Squared Error Decomposition
confidence 95% · we derive a mean-squared-error decomposition that separates model bias, estimation bias, model variance, and estimation variance for PD and ALE.
Cross-validation → reduces → Model Variance
confidence 90% · Cross-validation-based estimation is a promising approach that reduces the model variance component
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Global feature effects such as partial dependence (PD) and accumulated local effects (ALE) plots are widely used to interpret black-box models. However, they are only estimates of true underlying effects, and their reliability depends on multiple sources of error. Despite the popularity of global feature effects, these error sources are largely unexplored. In particular, the practically relevant question of whether to use training or holdout data to estimate feature effects remains unanswered. We address this gap by providing a systematic, estimator-level analysis that disentangles sources of bias and variance for PD and ALE. To this end, we derive a mean-squared-error decomposition that separates model bias, estimation bias, model variance, and estimation variance, and analyze their dependence on model characteristics, data selection, and sample size. We validate our theoretical findings through an extensive simulation study across multiple data-generating processes, learners, estimation strategies (training data, validation data, and cross-validation), and sample sizes. Our results reveal that, while using holdout data is theoretically the cleanest, potential biases arising from the training data are empirically negligible and dominated by the impact of the usually higher sample size. The estimation variance depends on both the presence of interactions and the sample size, with ALE being particularly sensitive to the latter. Cross-validation-based estimation is a promising approach that reduces the model variance component, particularly for overfitting models. Our analysis provides a principled explanation of the sources of error in feature effect estimates and offers concrete guidance on choosing estimation strategies when interpreting machine learning models.
Tags
Links
- Source: https://arxiv.org/abs/2603.15057v2
- Canonical: https://arxiv.org/abs/2603.15057v2
Trouble viewing inline? Open PDF directly →
Full Text
84,732 characters extracted from source content.
Expand or collapse full text
11institutetext: LMU Munich, Munich, Germany 11email: timo.heiss,giuseppe.casalicchio@stat.uni-muenchen.de 22institutetext: Munich Center for Machine Learning (MCML) Analyzing Error Sources in Global Feature Effect Estimation Timo HeiSS Coco Bögel Bernd Bischl Giuseppe Casalicchio Abstract Global feature effects such as partial dependence (PD) and accumulated local effects (ALE) plots are widely used to interpret black-box models. However, they are only estimates of true underlying effects, and their reliability depends on multiple sources of error. Despite the popularity of global feature effects, these error sources are largely unexplored. In particular, the practically relevant question of whether to use training or holdout data to estimate feature effects remains unanswered. We address this gap by providing a systematic, estimator-level analysis that disentangles sources of bias and variance for PD and ALE. To this end, we derive a mean-squared-error decomposition that separates model bias, estimation bias, model variance, and estimation variance, and analyze their dependence on model characteristics, data selection, and sample size. We validate our theoretical findings through an extensive simulation study across multiple data-generating processes, learners, estimation strategies (training data, validation data, and cross-validation), and sample sizes. Our results reveal that, while using holdout data is theoretically the cleanest, potential biases arising from the training data are empirically negligible and dominated by the impact of the usually higher sample size. The estimation variance depends on both the presence of interactions and the sample size, with ALE being particularly sensitive to the latter. Cross-validation-based estimation is a promising approach that reduces the model variance component, particularly for overfitting models. Our analysis provides a principled explanation of the sources of error in feature effect estimates and offers concrete guidance on choosing estimation strategies when interpreting machine learning models. 1 Introduction Many machine learning models are black boxes whose internal structure is not, or only partly, compatible with human reasoning, complicating explanations of both individual predictions and overall model behavior. This lack of transparency is especially problematic in high-stakes domains such as healthcare, law, and finance, where decisions must be transparent. To address this challenge, the field of eXplainable AI (XAI) has proposed a wide range of methods to explain machine learning models [molnar_interpretable_2022]. However, these methods must be used correctly to avoid misleading conclusions, as there are many pitfalls to be aware of [molnar_general_2022]. Global feature effect methods such as partial dependence (PD) [friedman_greedy_2001] and accumulated local effects (ALE) [apley_visualizing_2020] visualize how one or more features affect predictions. In practice, they are estimated from finite data, and their reliability depends on various error sources. Despite their widespread adoption, the error components of feature effect estimates remain largely unexplored. Prior works focus on extrapolation under feature dependence [apley_visualizing_2020, gkolemis_rhale_2023], aggregation bias [gkolemis_rhale_2023, goldstein_peeking_2015, herbinger_repid_2022, herbinger_decomposing_2024], or quantifying uncertainty [apley_visualizing_2020, cook_explaining_2024, moosbauer_explaining_2021]. A formal bias–variance decomposition w.r.t. a true underlying effect has only been derived for the theoretical PD [molnar_relating_2023], leaving estimator-level errors introduced by finite data largely unaddressed. A related practical question is whether to estimate explanations using training or holdout data. This has been studied for methods like permutation feature importance (PFI) [molnar_general_2022], mean decrease in impurity (MDI), and SHAP [loecher_debiasing_2022], but remains open for PD and ALE. While most works compute feature effects on training data [apley_visualizing_2020, friedman_greedy_2001, greenwell_pdp_2017, molnar_interpretable_2022], other works use holdout data [molnar_relating_2023]. Practitioners still debate whether to estimate PD and ALE on training or holdout data (see §2), trading larger training sample sizes against potential overfitting bias. These open questions point to a lack of an estimator-level understanding of error in global feature effect estimation, which our work addresses. Our main contributions are: • We provide the first estimator-level analysis of PD and ALE, deriving a full mean squared error (MSE) decomposition that separates model bias, estimation bias, model variance, and estimation variance. • We theoretically analyze these components, showing how sample size and interactions affect estimation bias and variance differently for PD and ALE, and formally relating remaining bias and variance to model bias and variance. • We empirically validate our findings in an extensive simulation study across multiple data-generating processes, learners, sample sizes, and estimation strategies (training, validation, and cross-validation (CV)), using dedicated estimators for the error components. We find negligible bias differences between training and holdout data, a strong sample-size effect – especially for ALE – and that CV is often preferable due to variance reduction. 2 Related Work Many issues in feature effects have been analyzed. The PD plot [friedman_greedy_2001] is known to suffer from extrapolation under dependent features by evaluating the model on implausible feature combinations [molnar_general_2022]. ALE plots avoid exactly that issue [apley_visualizing_2020]. Another issue is aggregation bias: global effects can obscure interaction-induced heterogeneity. Individual conditional expectation (ICE) curves [goldstein_peeking_2015] and RHALE [gkolemis_rhale_2023] visualize this heterogeneity, while regional methods like REPID [herbinger_repid_2022] and GADGET [herbinger_decomposing_2024] report effects in regions with reduced heterogeneity. Several works quantify uncertainty via variance of feature effects with model-specific PD confidence bands existing for probabilistic models [moosbauer_explaining_2021], as well as model-agnostic approaches that consider PD or ALE across multiple model fits [cook_explaining_2024, molnar_relating_2023, apley_visualizing_2020]. Only a few works explicitly study bias–variance trade-offs in feature effect estimators. RHALE [gkolemis_rhale_2023] optimizes ALE binning to balance bias and variance but does not relate this trade-off to a ground-truth effect. In [gkolemis_dale_2023], bias and variance of the proposed DALE estimator w.r.t. ALE are analyzed. In [chang_accelerated_2025], variance reduction and consistency of their proposed ALE estimator A2D2E are shown, and the error of PD, ALE, and A2D2E estimators against a derivative-based ground-truth are compared in simulations. Moreover, [molnar_relating_2023] formalize PD as a statistical estimator of a target estimand, derive a formal MSE decomposition for the theoretical PD, and propose variance estimators for confidence intervals. In contrast, we provide the first full estimator-level MSE decomposition of empirical PD and ALE w.r.t. their corresponding ground-truth feature effects, and analyze the error and its components both theoretically and empirically. For many interpretability methods, computing explanations on training vs. holdout data can affect conclusions: it matters for loss-based methods like permutation feature importance (PFI) [molnar_general_2022] and for mean decrease in impurity (MDI) and SHAP, which can be biased on training data [loecher_debiasing_2022]. For feature effects, this issue is largely unstudied. PD [friedman_greedy_2001], ALE [apley_visualizing_2020], and common references and software often use training data without justification [greenwell_pdp_2017, molnar_interpretable_2022], while others use holdout data [molnar_relating_2023]. Practitioners likewise disagree: some prefer training data for more reliable estimates due to larger sample size111https://forums.fast.ai/t/partial-dependence-plot/98465 (last accessed: 01/28/2026), others prefer holdout data as “the model may overfit”222https://github.com/sosuneko/PDPbox/issues/68 (01/28/2026), and some use both to diagnose distribution differences.333https://w.mathworks.com/help/stats/use-partial-dependence-plots-to-interpret-regression-models-trained-in-regression-learner-app.html (01/28/2026) To date, no systematic study addresses this. We do so by comparing different estimation strategies (training/holdout/CV) through bias-variance decomposition. 3 Notation & Background Assume a data-generating process that is characterized by a joint distribution ℙXYP_XY over features X=(X1,…,Xp)⊤X=(X_1,…,X_p) and target Y, with true underlying function f(x)=[Y∣X=x]f(x)=E[Y X=x]. A random dataset D=(X(i),Y(i))i=1nD=\(X^(i),Y^(i))\_i=1^n consists of n i.i.d. samples from ℙXYP_XY. Concrete realizations =((i),y(i))i=1nD=\(x^(i),y^(i))\_i=1^n are the training, test, and validation sets trainD_train, testD_test, and valD_val. A learning algorithm induces a fitted model f f on trainD_train, viewed either as a fixed function or as a random variable with distribution ℙFP_F due to training-sample and algorithmic randomness. Model performance is measured by the risk RL(f^)=XY[L(Y,f^(X))]R_L( f)=E_XY[L(Y, f(X))] with a point-wise loss L, and can be estimated on an n-sized dataset D as empirical risk Remp(f^;)=1n∑i=1nL(y(i),f^((i))).R_emp( f;D)= 1n _i=1^nL (y^(i), f(x^(i)) ). For trainD_train, this measures in-sample, for testD_test or valD_val out-of-sample error. Expectations and variances w.r.t. f^∼ℙF f _F are denoted by subscript F, D∼ℙXYD _XY by D, and (X,Y)∼ℙXY(X,Y) _XY by XYXY (analogously for marginal and conditionals). Note that D and f f are independent (denoted D⟂f^D f) only if D is not used for training, i.e., only for holdout data. We consider feature effects for a subset of features XSX_S with S⊆1,…,pS \1,…,p\. Throughout this work, we restrict attention to a single feature of interest (|S|=1|S|=1) and denote the complement feature set by XS¯X_ S. Definition 1(PD [friedman_greedy_2001]) Let h:→h:X be a prediction function and let XSX_S denote a feature of interest. The PD of h w.r.t. XSX_S is defined as PDh,S(xS):=XS¯[h(xS,XS¯)]=∫S¯h(xS,S¯)ℙXS¯(S¯).PD_h,S(x_S):=E_X_ S\! [h(x_S,X_ S) ]= _X_ Sh(x_S,x_ S)\,dP_X_ S(x_ S). (1) Given an n-sized dataset D, it can be estimated via Monte Carlo integration as PD^h,S(xS):=1n∑i=1nh(xS,S¯(i))=1n∑i=1nhS(i)(xS) PD_h,S(x_S):= 1n _i=1^nh(x_S,x_ S^(i))= 1n _i=1^nh^(i)_S(x_S). Here, hS(i)(xS)=h(xS,S¯(i))h^(i)_S(x_S)=h(x_S,x_ S^(i)) are the ICE curves [goldstein_peeking_2015]. PDf^,SPD_ f,S denotes the theoretical PD of the model f f, and PD^f^,S PD_ f,S is its estimator. PDf,SPD_f,S is the ground-truth PD. In practice, PD is visualized using a grid of G feature values xS(g)g=1G\x_S^(g)\_g=1^G. Quantile-based grids rather than equidistant ones are recommended in [molnar_general_2022]. Definition 2(ALE [apley_visualizing_2020]) Let h:→h:X be a prediction function and let XSX_S denote a feature of interest. The uncentered ALE of h w.r.t. XSX_S is defined as ALE~h,S(xS)=limK→∞∑k=1kSK(xS)XS¯|XS∈ISK(k)[Δh,SK(k,XS¯)], ALE_h,S(x_S)= _K→∞ _k=1^k_S^K(x_S)E_X_ S\,|\,X_S∈ I_S^K(k) [ _h,S^K(k,X_ S) ], (2) where feature XSX_S is partitioned into K intervals ISK(k)=(zk−1,SK,zk,SK]k=1K\I_S^K(k)=(z_k-1,S^K,z_k,S^K]\_k=1^K, kSK(xS)k_S^K(x_S) denotes the index of the interval into which a value xSx_S falls, and the maximum interval width converges to zero as K→∞K→∞. The finite differences in the k-th interval are given by Δh,SK(k,S¯):=h(zk,SK,S¯)−h(zk−1,SK,S¯) _h,S^K(k,x_ S):=h(z_k,S^K,x_ S)-h(z_k-1,S^K,x_ S). The uncentered ALE for a finite dataset D and finite K can be estimated as follows, where nSK(k)n_S^K(k) denotes the number of observations in the k-th interval: ALE~^h,S(xS)=∑k=1kSK(xS)1nSK(k)∑i:xS(i)∈ISK(k)[Δh,SK(k,S¯(i))] ALE_h,S(x_S)= _k=1^k_S^K(x_S) 1n_S^K(k) _i:x_S^(i)∈ I_S^K(k) [ _h,S^K(k,x_ S^(i)) ]. Centered versions ALEh,SALE_h,S and ALE^h,S ALE_h,S are obtained by subtracting a constant such that they have zero mean w.r.t. the marginal distribution of XSX_S or the empirical distribution of xS(i)i=1n\x_S^(i)\_i=1^n, respectively. For the grid zk,SKk=0K\z_k,S^K\_k=0^K, empirical quantiles of xS(i)i=1n\x_S^(i)\_i=1^n are recommended [apley_visualizing_2020]. For readability, we omit the superscript K when K is fixed, and suppress the subscript S for PD and ALE. In [molnar_relating_2023], it is shown that the MSE of the theoretical PD of f f w.r.t. to the theoretical ground-truth PD can be decomposed into squared bias and variance: F[(PDf(xS)−PDf^(xS))2]=(PDf(xS)−F[PDf^(xS)])2+VarF[PDf^(xS)].E_F[(PD_f(x_S)-PD_ f(x_S))^2]=(PD_f(x_S)-E_F[PD_ f(x_S)])^2+Var_F[PD_ f(x_S)]. (3) The bias term relates to systematic model bias, and the variance term captures variability across model fits. For the empirical estimator PD^f PD_ f, Molnar et al. [molnar_relating_2023] further argue that Monte Carlo integration introduces an additional source of variance. Moreover, when estimated on holdout data, PD^f PD_ f is unbiased w.r.t. PDf^PD_ f, and unbiasedness of the model implies unbiasedness of the PD.444Proofs can be found in Appendices C & D of the arXiv version of [molnar_relating_2023]. While an analogous decomposition for ALE is not available, consistency results for the ALE estimator exist in [apley_visualizing_2020]. To this end, they define a population version of the binned uncentered ALE for a prediction function h as ALE~hK(xS):=∑k=1kS(xS)XS¯|XS∈IS(k)[Δh(k,XS¯)] ALE_h^K(x_S):= _k=1^k_S(x_S)E_X_ S\,|\,X_S∈ I_S(k)\! [ _h(k,X_ S) ]. For the same K fixed bins, ALE~^h ALE_h on an i.i.d. D∼ℙXYD _XY converges to ALE~hK ALE_h^K pointwise as n→∞n→∞ almost surely, under mild integrability conditions. Moreover, ALE~hK ALE_h^K converges pointwise to ALE~h ALE_h as the bin resolution K→∞K→∞. Thus, ALE~^h ALE_h is jointly consistent for ALE~h ALE_h if n grows sufficiently fast relative to K.555These consistency results can be found in Theorem 3 of the arXiv version of [apley_visualizing_2020]. 4 Theoretical Considerations & Estimators For our theoretical analysis, we adopt several assumptions listed in §0.A.1. 4.1 Full Error Decomposition of the PD Estimator Previous work [molnar_relating_2023] considers only the error decomposition of PDf^PD_ f. Since PDf^PD_ f cannot be determined when ℙXS¯P_X_ S is unknown and is estimated via Monte Carlo integration, we instead study the estimator’s error. While [molnar_relating_2023] notes that this adds a variance term, we formally derive the MSE decomposition at fixed xSx_S, integrating over both f^∼ℙF f _F and the dataset D∼ℙXYD _XY used for PD estimation: FD∣f^[(PDf(xS)−PD^f^(xS))2]=(PDf(xS)−FD∣f^[PD^f^(xS)])2+VarFD∣f^[PD^f^(xS)]+FVarD∣f^[PD^f^(xS)]. splitE_FE_D f [(PD_f(x_S)- PD_ f(x_S))^2 ]=(PD_f(x_S)-E_FE_D f[ PD_ f(x_S)])^2\\ +Var_FE_D f[ PD_ f(x_S)]+E_FVar_D f[ PD_ f(x_S)]. split (4) The proof is given in §0.A.2. The decomposition has three distinct terms, which we will analyze in more detail below: the squared bias and two variances. Bias. The bias of the PD estimator (cf. first term in Eq. (4)) decomposes as: FD∣f^[PD^f^(xS)]−PDf(xS)=F[D∣f^[PD^f^(xS)]−PDf^(xS)]+(F[PDf^(xS)]−PDf(xS)) splitE_FE_D f[ PD_ f(x_S)]-PD_f(x_S)=\\ E_F[E_D f[ PD_ f(x_S)]-PD_ f(x_S)]+(E_F[PD_ f(x_S)]-PD_f(x_S)) split (5) by adding and subtracting F[PDf^,S(xS)]E_F[PD_ f,S(x_S)]. The first part vanishes by the unbiasedness of the PD estimator w.r.t. the theoretical model PD [molnar_relating_2023], if only the data used for Monte Carlo integration in the PD estimation is independent of the model f f. This is true for holdout data, but not necessarily for training data. For the second part, exchanging expectations (Fubini) yields (proof in §0.A.3): F[PDf^(xS)]−PDf(xS)=XS¯[F[f^(xS,XS¯)]−f(xS,XS¯)].E_F\! [PD_ f(x_S) ]-PD_f(x_S)=E_X_ S\! [E_F\! [ f(x_S,X_ S) ]-f(x_S,X_ S) ]. (6) Consequently, the term reduces to the model’s bias averaged over ℙXS¯P_X_ S. Thus, for estimation on holdout data, the PD estimator’s bias reduces to the average model bias, whereas estimation on training data may introduce additional bias. Model variance. The second term in Eq. (4) reflects the variance of the PD estimator w.r.t. the model distribution ℙFP_F. When the PD is estimated on holdout data, it reduces to VarF[D∣f^[PD^f^(xS)]]=VarF[PDf^(xS)]Var_F[E_D f[ PD_ f(x_S)]]=Var_F[PD_ f(x_S)] by the unbiasedness of the PD estimator w.r.t. the theoretical model PD. This is exactly the variance of the theoretical PD in Eq. (3). By exchanging expectations (Fubini/Tonelli) and applying Jensen’s inequality, we obtain an upper bound (proof in §0.A.4): VarF[PDf^(xS)]≤XS¯VarF[f^(xS,XS¯)].Var_F[PD_ f(x_S)] _X_ S\!Var_F[ f(x_S,X_ S)]. (7) Thus, the theoretical PD variance at xSx_S is controlled by the average (pointwise) model variance w.r.t. the marginal distribution of XS¯X_ S. Estimation variance. The third term in Eq. (4) captures the variance w.r.t. the samples used to estimate the PD via Monte Carlo integration, which equals: FVarD∣f^[PD^f^(xS)]=1nFVarD∣f^[f^(xS,XS¯)]E_FVar_D f[ PD_ f(x_S)]= 1nE_FVar_D f[ f(x_S,X_ S)] (8) for any XS¯∼D|f^X_ S D| f. The proof is given in §0.A.5. Consequently, the pointwise estimation variance of the PD depends on the sample size n and decreases at the rate (1/n)O(1/n). It also depends on the expected variance of the ICE curves at xSx_S. Centering the PD via PD^f^(xS)−D∣f^[PD^f^(XS)] PD_ f(x_S)-E_D f[ PD_ f(X_S)] yields the variance of the centered ICE curves VarD∣f^[f^(xS,XS¯)−D∣f^[f^(XS,XS¯)]]Var_D f[ f(x_S,X_ S)-E_D f[ f(X_S,X_ S)]] instead (proof analogous to §0.A.5). Since centered ICE curves fulfill local decomposability [herbinger_decomposing_2024], their variance at xSx_S is solely due to interactions involving feature XSX_S. Thus, the estimation variance of the centered PD is zero when XSX_S has no interactions in f f. An overview of all four derived error components is provided in Fig. 1. Figure 1: Conceptual overview of the four error components. 4.2 Full Error Decomposition of the ALE Estimator We now provide an error analysis for ALE, which is missing in the current literature. The MSE of ALE^f ALE_ f w.r.t. ALEfALE_f can be decomposed into bias and variance analogous to PD in Eq. (4) (proof in §0.A.2, analogous to the one for PD). The only property required for the argument is that all expectations exist and that f, and hence ALEfALE_f, is non-random w.r.t. to ℙFP_F. For the theoretical analysis of the error components, we focus on uncentered ALE, as centering is linear post-processing that only affects the offset. All sources of bias and variance originate in the uncentered ALE and propagate deterministically under centering. Bias. Similar to PD, we can decompose the bias of ALE further into: FD∣f^[ALE~^f^(xS)]−ALE~f(xS)=F[D∣f^[ALE~^f^(xS)]−ALE~f^(xS)]+(F[ALE~f^(xS)]−ALE~f(xS)), splitE_FE_D f[ ALE_ f(x_S)]- ALE_f(x_S)=\\ E_F\! [E_D f[ ALE_ f(x_S)]- ALE_ f(x_S) ]+ (E_F[ ALE_ f(x_S)]- ALE_f(x_S) ), split (9) by adding and subtracting F[ALE~f^(xS)]E_F[ ALE_ f(x_S)]. The first part can be viewed as the average “estimation bias” and can be further decomposed into: F[D∣f^[ALE~^f^(xS)]−ALE~f^(xS)]=F[D∣f^[ALE~^f^(xS)]−ALE~f^K(xS)]+F[ALE~f^K(xS)−ALE~f^(xS)]. splitE_F[E_D f[ ALE_ f(x_S)]- ALE_ f(x_S)]=\\ E_F[E_D f[ ALE_ f(x_S)]- ALE_ f^K(x_S)]+E_F[ ALE_ f^K(x_S)- ALE_ f(x_S)]. split (10) The first term is due to estimation on finite samples. For ALE estimation on holdout data, under mild regularity (Assumptions (i) and (v)) and nS(k)>0n_S(k)>0 for all relevant bins, the ALE estimator is unbiased w.r.t. the binned population ALE and the term becomes zero (proof in §0.A.6). To the second term, we refer as “discretization bias”: the term inside FE_F goes to 0 as K→∞K→∞ (pointwise in f f) by the definition of ALE (Eq. (2)). Convergence of the entire expectation term can be shown by applying the dominated convergence theorem exactly as in §0.A.7. For the second part in Eq. (9), we show in §0.A.7 that it arises from an (infinite) sum of the local average biases in the finite differences of f f w.r.t. those of f. Thus, it vanishes when the model’s conditional finite differences are unbiased w.r.t. those of f. In particular, this is satisfied if f f is unbiased w.r.t. f for all x. Model variance. When ALE is estimated on holdout data, it holds that VarF[D∣f^[ALE~^f^(xS)]]=VarF[ALE~f^K(xS)]Var_F[E_D f[ ALE_ f(x_S)]]=Var_F[ ALE_ f^K(x_S)] by the ALE estimator’s unbiasedness (see above). For this expression, we obtain an upper bound for the variance of ALE w.r.t. the model distribution (proof in §0.A.8): VarF[ALE~f^K(xS)]≤kS(xS)∑k=1kS(xS)XS¯∣XS∈IS(k)VarF[Δf^(k,XS¯)].Var_F\! [ ALE_ f^K(x_S) ]≤ k_S(x_S)\!\! _k=1^k_S(x_S)E_X_ S X_S∈ I_S(k)Var_F\! [ _ f(k,X_ S) ]. (11) Thus, the theoretical ALE variance at xSx_S is controlled by the local average variability of finite differences across models along the intervals up to xSx_S. Estimation variance. By conditioning on the samples’ bin assignments, B=(B(1),…,B(n))B=(B^(1),…,B^(n)) with B(i)=kS(XS(i))B^(i)=k_S(X_S^(i)), and using the law of total variance, the estimation variance of ALE, F[VarD∣f^[ALE~^f^(xS)]]E_F[Var_D f[ ALE_ f(x_S)]], decomposes into F[VarB∣f^[D∣f^,B[ALE~^f^(xS)]]]+F[B∣f^[VarD∣f^,B[ALE~^f^(xS)]]]. _F\! [Var_B f\! [E_D f,B\! [ ALE_ f(x_S) ] ] ]+E_F\! [E_B f\! [Var_D f,B\! [ ALE_ f(x_S) ] ] ]. The first term captures variance from random bin assignments along XSX_S, i.e., variability in the expected ALE estimate due to random sample allocation to bins. For the second term, assuming nS(k)>0n_S(k)>0 for all relevant bins k≤kS(xS)k≤ k_S(x_S): FB∣f^VarD∣f^,B[ALE~^f^(xS)]=∑k≤kS(xS)FB∣f^[1nS(k)σk2(f^)]E_FE_B fVar_D f,B\! [ ALE_ f(x_S) ]= _k≤ k_S(x_S)E_FE_B f [ 1n_S(k) _k^2( f) ] (12) with σk2(f^):=VarXS¯∣XS∈IS(k),f^[Δf^(k,XS¯)] _k^2( f):=Var_X_ S X_S∈ I_S(k), f[ _ f(k,X_ S)]. The proof is given in §0.A.9. Thus, it scales with the expected inverse number of observations per bin. Assuming deterministic equal-frequency binning for simplicity, this factor becomes Kn Kn. It also depends on the local variance in the estimated finite differences. Finite differences also fulfill local decomposability [herbinger_decomposing_2024]. Thus, the variance σk2(f^) _k^2( f) is only due to interaction effects involving feature XSX_S and is zero when XSX_S has no interaction effects in f f. 4.3 Estimators of the Error Components To investigate the error components empirically in §5-6, we propose estimators for each. FE_F and VarFVar_F can be estimated by averaging over multiple models f^(m) f^(m) of the same learner fitted to M different training sets, each independently sampled from ℙXYP_XY. Translating the standard estimators for MSE, bias, and variance [morris_using_2019] to our setting yields the estimators in Tab. 1. The variance estimator is the same as in [molnar_relating_2023], capturing the variance in both model fits and Monte Carlo integration. To estimate the estimation variance separately, consider R Monte Carlo iterations rD_r to estimate the feature effect PD^[r] PD^[r], again for multiple model refits. While Tab. 1 is based on PD, estimators for ALE follow analogously. Table 1: Estimators of the feature effect error components at xSx_S. Component Estimator MSE MSE^(xS)=1M∑m=1M(PDf(xS)−PD^f^(m)(xS))2 MSE(x_S)= 1M _m=1^M (PD_f(x_S)- PD_ f^(m)(x_S) )^2 Bias Bias^(xS)=PDf(xS)−1M∑m=1MPD^f^(m)(xS) Bias(x_S)=PD_f(x_S)- 1M _m=1^M PD_ f^(m)(x_S) Variance Var^(xS)=1M−1∑m=1M(PD^f^(m)(xS)−1M∑m′=1MPD^f^(m′)(xS))2 Var(x_S)= 1M-1 _m=1^M ( PD_ f^(m)(x_S)- 1M _m =1^M PD_ f^(m )(x_S) )^2 Estimation variance Var^Est(xS)=1M(R−1)∑m=1M∑r=1R(PD^f^(m)[r](xS)−1R∑r′=1RPD^f^(m)[r′](xS))2 Var_Est(x_S)= 1M(R-1) _m=1^M _r=1^R ( PD_ f^(m)^[r](x_S)- 1R _r =1^R PD_ f^(m)^[r ](x_S) )^2 These estimators are unbiased according to standard statistical results. The Monte Carlo standard errors of these estimates can again be estimated as a function of the number of repetitions [morris_using_2019]. Importantly, most estimators in Tab. 1 are impractical as they require knowledge of the ground-truth. We will rather use them to understand the behavior of feature effect estimation. 5 Experimental Set-Up We now empirically validate our findings from §4 with an extensive simulation study. While the theoretical analysis requires holdout data at multiple points (e.g., to avoid estimation bias), we aim to gain deeper insights into what empirically happens when this is violated. For this, we compare feature effect errors across different estimation strategies: on training data, validation data, and via CV. We break this overarching goal down into three specific research questions: RQ1: How does overfitting empirically affect MSE, bias, and variance of PD and ALE when estimated on training vs. validation data vs. CV? RQ2: How do the variance sources (model and estimation variance) empirically behave for PD and ALE on training vs. validation data vs. in CV? RQ3: How does sample size affect estimation error for PD and ALE empirically? Data settings. We consider three settings of varying complexity. Ground-truth functions and feature structures are given in Tab. 2. All settings include two independent dummy features. The first setting has correlations and interactions, the second different non-linearities, and the third real-world relevance.666This dataset is based on the physics-grounded Feynman equation I.29.16 for wave interference, addressing the limited realism of standard test functions. Features are generated i.i.d., and the target as y=f()+εy=f(x)+ with i.i.d. ε∼(0,σε2) (0, _ ^2) and a signal-to-noise ratio of 55.777The scale parameter is set such that σ^Y/σε=5 σ_Y/ _ =5. We consider two sample sizes n=1250n=1250 and n=10000n=10000. Table 2: Data settings with underlying functions and feature structure. Setting Function Correlation Simple-Normal-Correlated f1()=x1+x222+x1x2 array[]lf_1(x)=x_1+ x_2^22+x_1x_2 array X1,…,X4∼N(0,1),ρ12=0.9,ρij=0∀(i,j)≠(1,2),i≠j array[]lX_1,…,X_4 N (0,1 ), _12=0.9,\\ _ij=0\;∀(i,j)≠(1,2),\ i≠ j array Friedman1 [friedman_multivariate_1991] f2()=10sin(πx1x2)+20(x3−12)2+10x4+5x5 array[]lf_2(x)=10 (π x_1x_2)\\ +20(x_3- 12)^2+10x_4+5x_5 array X1,…,X7∼i.i.d.U(0,1)X_1,…,X_7 i.i.d. U(0,1) Feynman I.29.16 [matsubara_rethinking_2024] f3(,)=x12+x22+2x1x2cos(θ1−θ2) array[]lf_3(x, θ)= x_1^2+x_2^2+2x_1x_2 ( _1- _2) array X1,X2∼LogU(0.1,10),θ1,θ2∼U(0,2π),D1,D2∼U(0,1), all indep. array[]lX_1,X_2 (0.1,10),\\ _1, _2 (0,2π),\\ D_1,D_2 (0,1), all indep. array Models. We consider a generalized additive model (GAM) with spline bases for main and pairwise interaction effects, and XGBoost as learners. Each of them is once configured with “optimally tuned” (OT) hyperparameters and once with hyperparameters chosen to overfit (OF). These hyperparameters are pre-selected per setting and sample size, i.e., carefully hand-picked for OF (e.g., small penalty / large learning rate) and tuned on separate data samples for OT. For details on the hyperparameters and model performances, see §0.B.1. Feature effect estimation. We estimate the feature effects PD^f PD_ f and ALE^f ALE_ f per feature and model by (a) training a model on all n samples and estimating the feature effect on the same n samples, (b) splitting the n samples into 80%80\% train and 20%20\% validation set, fitting a model on the training set and estimating the effect on the validation set, and (c) a CV-based estimation strategy on all n samples.888We use 5-fold CV: in each iteration, a model is fitted on four folds, effects are estimated on the held-out fold, and the five resulting effects are averaged pointwise. We compute feature effects all at the same 100 grid points, defined by the theoretical quantiles of ℙXSP_X_S for comparability. Additionally, we center the curves after estimation, and omit the first and last grid point to avoid boundary effects, particularly for ALE. To enable ground-truth comparisons, we additionally construct ground-truth effect estimators in the same manner, but on 10,00010,000 fresh samples and directly on f. This eliminates the discretization bias from our error, as the ground-truth and estimate use the same intervals.999Except for discretization bias, this estimate is unbiased w.r.t. the true theoretical effect (cf. §4), and adds negligible estimation variance at this sample size (cf. §6.3). Experiments. To address RQ1, we compute the MSE, bias, and variance of all estimated model feature effects w.r.t. the estimated ground-truth effects via the estimators in Tab. 1 with PD^f PD_f for PDfPD_f (analog. for ALE). We repeat each setting-size-model combination M=30M=30 times to estimate the error terms. Each repetition involves drawing new data, fitting the models, and estimating the effects. We report MSE, bias, and variance averaged over the grid points. To address RQ2, we estimate the estimation variance according to Tab. 1. In each iteration m∈1,…,Mm∈\1,…,M\, we fix the trained models, draw R=30R=30 new data sets, and estimate the effects with them for the fixed models. By subtracting this from the total variance (cf. RQ1), we estimate the model variance. For this analysis, we focus on XGBoost (OF & OT) and a sample size of n=1250n=1250. To address RQ3, we compute the MSE between analytical ground-truth effects (available for the first two settings) and estimated ground-truth effects across 50 different sample sizes ranging from 10110^1 to 10610^6 (on log scale) with 50 repetitions per size. This isolates the estimation error, as no model is involved. Reproducibility. Experiments are implemented in Python and use fixed random seeds. We release all experimental code and raw results on GitHub (link). 6 Empirical Results We report results for a single setting and a representative feature subset. Further results are provided in §0.B.2 and consistent with the findings reported here. 6.1 Bias-Variance-Analysis PD decomposition. Our empirical results on SimpleNormalCorrelated in Tab. 3 reveal that the MSE of the PD estimator is lowest mostly for CV- and training-set-based estimation. Bias is mostly similar across the estimation strategies, and we observe no systematic trends w.r.t. estimation strategy or sample size. Thus, a potential bias introduced by estimation on training data that could not be ruled out in our theoretical analysis (§4.1) appears empirically negligible. Variance is generally lowest for CV-based estimation. We hypothesize that this is due to two effects: (1) model variance may decrease as CV averages out fitting variability across multiple models, and (2) estimation variance may decrease compared to estimation on a single validation set due to increased effective sample size. Variance is generally slightly higher for validation than for training-set-based estimation, likely due to the smaller sample size. The higher variance of overfitting models reflects in higher PD variance. These empirical results agree with our theoretical findings. An additional finding is that, for well-generalizing models, differences in MSE across estimation strategies are mostly negligible, while for overfitting models, CV-based estimation yields a substantial reduction. Table 3: Results for PD on SimpleNormalCorrelated averaged over 100 grid points. Bold numbers are minimum per metric-feature-model-size-combination. Feature x_1 x_2 x_3 Metric MSE Bias Var MSE Bias Var MSE Bias Var n=1250n=1250 GAM_OF train 0.0695 0.0296 0.0710 0.0699 0.0276 0.0716 0.0090 0.0161 0.0091 val 0.0928 0.0206 0.0955 0.0849 0.0345 0.0865 0.0114 0.0190 0.0114 CV 0.0609 0.0264 0.0622 0.0562 0.0292 0.0573 0.0086 0.0157 0.0086 GAM_OT train 0.0051 0.0323 0.0042 0.0039 0.0247 0.0034 0.0005 0.0022 0.0005 val 0.0089 0.0336 0.0081 0.0093 0.0305 0.0086 0.0005 0.0040 0.0005 CV 0.0053 0.0356 0.0042 0.0041 0.0287 0.0034 0.0005 0.0022 0.0005 XGB_OF train 0.1780 0.3147 0.0817 0.4217 0.4505 0.2263 0.0010 0.0043 0.0010 val 0.1880 0.3242 0.0857 0.3888 0.3941 0.2416 0.0015 0.0098 0.0014 CV 0.1458 0.3114 0.0505 0.3043 0.4236 0.1291 0.0007 0.0042 0.0007 XGB_OT train 0.2807 0.5194 0.0113 0.1690 0.3941 0.0141 0.0014 0.0052 0.0014 val 0.3008 0.5375 0.0123 0.1666 0.3894 0.0155 0.0019 0.0073 0.0019 CV 0.2950 0.5351 0.0090 0.1691 0.3987 0.0104 0.0013 0.0051 0.0013 n=10000n=10000 GAM_OF train 0.1878 0.0674 0.1895 0.1995 0.0890 0.1982 0.0008 0.0050 0.0008 val 0.2248 0.0601 0.2289 0.2371 0.0678 0.2406 0.0011 0.0055 0.0011 CV 0.1777 0.0778 0.1775 0.1875 0.0996 0.1837 0.0008 0.0050 0.0008 GAM_OT train 0.0011 0.0206 0.0007 0.0012 0.0204 0.0008 0.0000 0.0009 0.0000 val 0.0020 0.0228 0.0016 0.0018 0.0222 0.0014 0.0001 0.0012 0.0001 CV 0.0012 0.0232 0.0007 0.0013 0.0230 0.0008 0.0000 0.0010 0.0000 XGB_OF train 0.0915 0.2638 0.0227 0.2625 0.4404 0.0710 0.0003 0.0029 0.0003 val 0.1029 0.2675 0.0324 0.2700 0.4465 0.0730 0.0004 0.0037 0.0004 CV 0.0828 0.2641 0.0135 0.2226 0.4277 0.0411 0.0002 0.0022 0.0002 XGB_OT train 0.1752 0.4142 0.0037 0.1210 0.3418 0.0043 0.0002 0.0022 0.0002 val 0.1802 0.4197 0.0042 0.1226 0.3439 0.0045 0.0003 0.0022 0.0003 CV 0.1806 0.4215 0.0030 0.1237 0.3471 0.0033 0.0002 0.0021 0.0002 ALE decomposition. Similar results for ALE in Tab. 4 show that the MSE is again lowest for training-set- or CV-based estimation. In contrast to PD, bias is often considerably lower for training data (largest set) when n is small, and generally decreases with increasing sample size, likely as the probability of nS(k)>0∀kn_S(k)>0\ ∀ k grows (required for unbiasedness, cf. §4.2). Again, the variance is mostly lowest for CV-based estimation, but now substantially higher on the smaller validation set, supporting our theoretical result that ALE (estimation) variance is more sensitive to sample size than PD. Generally, this confirms our theoretical findings and again shows that CV-based estimation is promising. Table 4: Results for ALE on SimpleNormalCorrelated averaged over 100 grid points. Bold numbers are minimum per metric-feature-model-size-combination. Feature x_1 x_2 x_3 Metric MSE Bias Var MSE Bias Var MSE Bias Var n=1250n=1250 GAM_OF train 0.0123 0.0205 0.0123 0.0108 0.0236 0.0106 0.0091 0.0161 0.0092 val 0.0277 0.0771 0.0225 0.0304 0.0711 0.0262 0.0129 0.0176 0.0131 CV 0.0178 0.0857 0.0108 0.0190 0.0858 0.0120 0.0077 0.0153 0.0077 GAM_OT train 0.0022 0.0100 0.0021 0.0019 0.0061 0.0019 0.0005 0.0022 0.0005 val 0.0148 0.0808 0.0086 0.0198 0.0918 0.0117 0.0005 0.0036 0.0005 CV 0.0103 0.0899 0.0023 0.0109 0.0890 0.0031 0.0004 0.0020 0.0004 XGB_OF train 0.1159 0.1343 0.1012 0.1120 0.1196 0.1011 0.0724 0.0279 0.0741 val 0.5238 0.1153 0.5281 0.6017 0.1415 0.6017 0.0848 0.0354 0.0864 CV 0.1531 0.0784 0.1520 0.1025 0.0803 0.0994 0.0118 0.0126 0.0121 XGB_OT train 0.0162 0.0875 0.0088 0.0174 0.0768 0.0119 0.0065 0.0077 0.0066 val 0.0471 0.1470 0.0263 0.0386 0.1094 0.0275 0.0033 0.0116 0.0033 CV 0.0315 0.1549 0.0077 0.0234 0.1158 0.0103 0.0020 0.0053 0.0020 n=10000n=10000 GAM_OF train 0.0011 0.0071 0.0011 0.0012 0.0071 0.0011 0.0008 0.0050 0.0008 val 0.0017 0.0075 0.0017 0.0017 0.0085 0.0017 0.0011 0.0057 0.0011 CV 0.0011 0.0068 0.0011 0.0011 0.0072 0.0011 0.0008 0.0050 0.0008 GAM_OT train 0.0002 0.0030 0.0002 0.0003 0.0045 0.0002 0.0000 0.0009 0.0000 val 0.0006 0.0037 0.0006 0.0005 0.0056 0.0005 0.0001 0.0012 0.0001 CV 0.0002 0.0029 0.0002 0.0003 0.0046 0.0002 0.0000 0.0010 0.0000 XGB_OF train 0.0164 0.0346 0.0157 0.0197 0.0520 0.0176 0.0103 0.0189 0.0103 val 0.0652 0.0416 0.0657 0.0780 0.0597 0.0770 0.0155 0.0166 0.0157 CV 0.0130 0.0180 0.0131 0.0147 0.0198 0.0148 0.0023 0.0098 0.0023 XGB_OT train 0.0046 0.0538 0.0017 0.0039 0.0401 0.0023 0.0005 0.0037 0.0005 val 0.0054 0.0544 0.0025 0.0047 0.0412 0.0031 0.0006 0.0037 0.0007 CV 0.0044 0.0555 0.0014 0.0034 0.0417 0.0017 0.0003 0.0027 0.0003 6.2 Variance Decomposition Analysis For variance decomposition into model and estimation variance, we consider the results on Friedman1 in Tab. 5, including a non-linear feature with interactions (X1X_1), a linear without (X4X_4), and a dummy feature (X7X_7). The estimation variance is constantly highest when feature effects are estimated on the smaller validation set. This is more pronounced for ALE, and we generally observe higher estimation variance than for PD. As hypothesized, CV (1) reduces model variance compared to the other estimation strategies, which is most pronounced for overfitting models (generally higher model variance). As expected, it also (2) reduces estimation variance compared to a single validation set. These empirical findings agree with our theoretical results. Although estimation variance is sometimes slightly higher for features with interactions (X1X_1), XGBoost may also learn interactions that are not in f, which is why this is not as clear as expected. Table 5: Decomposition results of total variance VarTotVar_Tot into model variance VarModVar_Mod and estimation variance VarEstVar_Est on Friedman1 averaged over 100 grid points. Bold numbers indicate estimation strategy with minimal variance. Feature x_1 x_4 x_7 Metric VarTotVar_Tot VarModVar_Mod VarEstVar_Est VarTotVar_Tot VarModVar_Mod VarEstVar_Est VarTotVar_Tot VarModVar_Mod VarEstVar_Est PD XGB_OF train 0.0862 0.0833 0.0029 0.1269 0.1240 0.0029 0.0027 0.0025 0.0002 val 0.1120 0.0974 0.0146 0.1811 0.1661 0.0150 0.0059 0.0046 0.0012 CV 0.0460 0.0432 0.0028 0.0814 0.0785 0.0029 0.0014 0.0011 0.0002 XGB_OT train 0.0217 0.0209 0.0008 0.0232 0.0231 0.0001 0.0071 0.0070 0.0000 val 0.0280 0.0241 0.0038 0.0283 0.0278 0.0006 0.0083 0.0082 0.0001 CV 0.0190 0.0182 0.0008 0.0205 0.0204 0.0001 0.0062 0.0062 0.0000 ALE XGB_OF train 1.0248 0.4286 0.5962 1.1178 0.4220 0.6958 0.1852 0.1556 0.0296 val 4.6600 1.5201 3.1399 4.9819 1.1399 3.8420 0.3004 0.1614 0.1390 CV 0.6573 0.0846 0.5727 0.8568 0.0858 0.7710 0.0262 -* 0.0287 XGB_OT train 0.0337 0.0229 0.0108 0.0340 0.0265 0.0075 0.0125 0.0113 0.0012 val 0.0937 0.0200 0.0737 0.1113 0.0463 0.0651 0.0154 0.0075 0.0080 CV 0.0299 0.0147 0.0151 0.0273 0.0143 0.0129 0.0067 0.0051 0.0016 * Due to instabilities in variance estimation for ALE, likely caused by unstable bin assignments, the estimated VarEstVar_Est sometimes exceeds the estimated VarTotVar_Tot. We omit VarModVar_Mod in these cases. 6.3 Effect of the Sample Size Our results for the effect of sample size are shown in Fig. 2 for Friedman1. For PD on holdout data, we know that there is no estimation bias (cf. §4.1) and Fig. 2(a) contains only estimation variance. We observe the expected polynomial decrease at roughly 1/n1/n for features with interactions (X1X_1) and a negligible error across sample sizes for features without interactions (X4X_4), as expected for centered PD. For ALE, we know from §4.2 that there are also estimation biases in addition to variance due to discretization and when nS(k)=0n_S(k)=0. With interactions (X1X_1), the observed estimation error in Fig. 2(b) is closer to K/nK/n for small sample sizes but gets closer to 1/n1/n as the sample size increases. For features without interactions, we observe a sharp drop in the estimation error at n=Kn=K, reducing it to the expected negligible level. This may be due to reduced estimation bias at this sample size, as each interval could, in principle, now contain at least one sample. (a) Mean estimation error for PD (b) Mean estimation error for ALE Figure 2: Mean estimation errors on Friedman1 for X1X_1 and X4X_4. For each sample size n, the variances are averaged over all grid points. Both axes are log-scale. 7 Conclusion & Future Work Summary. In this work, we presented an estimator-level analysis of global feature effect estimation for PD and ALE. We derived a full MSE decomposition that separates model bias, estimation bias, model variance, and estimation variance, and we analyzed these components theoretically. Our results show that the model bias component follows directly from systematic biases in the model f f for PD or in the finite differences of f f for ALE. For PD, the estimation bias is zero on holdout data. For ALE, it consists of a bias when nS(k)=0n_S(k)=0, and a discretization bias due to binning. For both PD and ALE, we derived upper bounds on the model variance in terms of pointwise variance of f f (PD) or its finite differences (ALE). Estimation variance is governed by (1) the sample size, scaling as 1/n1/n for PD and as the expected inverse bin counts for ALE, explaining ALE’s sensitivity in small-sample regimes, and (2) the variance of the ICE curves / finite differences. For centered PD and generally for ALE, the latter depends only on interactions with XSX_S. We validated our theoretical findings through simulations by comparing feature effect estimation on training and validation data and a CV-based strategy. Our empirical results showed that potential bias from estimating feature effects on the training data is negligible across our settings, including when models overfit (RQ1). Differences are instead dominated by sample-size effects: validation-set-based estimation yields higher variance (for ALE also higher bias), while training-set- and CV-based estimation attain the lowest MSE. Variance decomposition (RQ2) shows that the estimation variance is consistently highest on the smaller validation set and more pronounced for ALE, in line with our theoretical results. CV reduces model variance by averaging across fits (particularly beneficial for overfitting models) and estimation variance through a higher effective sample size compared to a single validation set. The observed estimation error (RQ3) confirms our theoretical results on sample size and interaction effects. Practical implications. While using holdout data is theoretically cleaner, our results indicate that training-set-based feature effect estimation is empirically safe and often preferable because it has a larger sample size. CV-based estimation emerges as a robust alternative, particularly for overfitting models. Although such models are typically identified via appropriate performance estimation and ruled out, this may not always be perfectly possible in all applications, in which cases CV-based feature effect estimation can be a safer option. Limitations & future work. Our empirical results are themselves estimates with errors from finite repetitions, which can lead to artifacts, such as estimation variance estimates exceeding total variance (cf. Tab. 5). Additionally, the empirical analysis considers only low-dimensional settings and two model classes. Further, our theoretical analysis leaves room for future research, including: (1) a tighter bias analysis for cases in which D⟂f^D f does not hold to formalize what happens for training-set-based estimation, (2) a formal theory for CV-based estimation results, (3) extensions of our analysis to distribution shifts between training and estimation data, which is another practically relevant issue. References Appendix Appendix 0.A Theoretical Evidence 0.A.1 Assumptions (i) Independence of feature effect definition and model. Feature effects (both PD and ALE) are defined using expectations independent of the randomness in the training of f^∼ℙF f _F, meaning XS¯⟂f^X_ S f for the theoretical PD / ALE. (i) Integrability of f f and f. For all features XSX_S and points xSx_S, the random variable f^(xS,XS¯) f(x_S,X_ S) is integrable and square-integrable under the joint law ℙ(F,XS¯)P_(F,X_ S), i.e., (F,XS¯)[|f^(xS,XS¯)|]<∞E_(F,X_ S)[| f(x_S,X_ S)|]<∞ and (F,XS¯)[f^(xS,XS¯)2]<∞E_(F,X_ S)[ f(x_S,X_ S)^2]<∞. In particular, this implies that for ℙFP_F-almost all realizations of f f we have XS¯∣f^[f^(xS,XS¯)2]<∞E_X_ S f[ f(x_S,X_ S)^2]<∞ by Tonelli’s theorem for conditional probabilities / stochastic kernels (since f^(xS,XS¯)2≥0 f(x_S,X_ S)^2≥ 0). Likewise XS¯∣f^[|f^(xS,XS¯)|]<∞E_X_ S f[| f(x_S,X_ S)|]<∞. Analogously, f(xS,XS¯)f(x_S,X_ S) is integrable and square-integrable under ℙXS¯P_X_ S, i.e., XS¯[|f(xS,XS¯)|]<∞E_X_ S[|f(x_S,X_ S)|]<∞ and XS¯[f(xS,XS¯)2]<∞E_X_ S[f(x_S,X_ S)^2]<∞. These conditions hold in particular for the special case that XS¯X_ S and f f are independent, e.g., under Assumption (i).They ensure that the relevant expectations and variances are well defined (i.p., PD, ALE, and their estimators) and justify interchanging the order of integration (Fubini/Tonelli) where needed. (i) Finite conditional probabilities. For all bins IS(k)I_S(k) considered, ℙ(XS∈IS(k))>0,P(X_S∈ I_S(k))>0, and the conditional law ℙXS¯∣XS∈IS(k)P_X_ S X_S∈ I_S(k) is well-defined. (iv) Uniform bounded total variation in xSx_S. For all features XSX_S, there exists a nonnegative random variable V(f^)V( f) such that, with probability 1 over f^∼ℙF f _F, esssupS¯TV(t↦f^(t,S¯))≤V(f^),F[V(f^)]<∞. *ess\,sup_x_ STV\! (t f(t,x_ S) )≤ V( f), _F[V( f)]<∞. (v) Integrability of finite differences. For all features XSX_S, all K and all k considered, the finite differences of f f are integrable and square-integrable w.r.t. the joint law given the k-th bin, i.e., under ℙ(F,X)∣XS∈ISK(k)P_(F,X) X_S∈ I_S^K(k). In other words, (F,X)∣XS∈ISK(k)[Δf^,SK(k,XS¯)2]<∞E_(F,X) X_S∈ I_S^K(k)[\, _ f,S^K(k,X_ S)^2]<∞, and analogously to (i) this implies XS¯∣XS∈ISK(k),f^[Δf^,SK(k,XS¯)2]<∞E_X_ S X_S∈ I_S^K(k), f[\, _ f,S^K(k,X_ S)^2]<∞ for ℙFP_F-almost all f f. For integrability, we likewise have (F,X)∣XS∈ISK(k)[|Δf^,SK(k,XS¯)|]<∞E_(F,X) X_S∈ I_S^K(k)[\,| _ f,S^K(k,X_ S)|]<∞. Note that this follows from (iv) plus the square-integrability condition in (i). As in (i), the same integrability conditions hold for Δf,SK(k,XS¯) _f,S^K(k,X_ S) (with f in place of f f). (vi) Existence and square-integrability of ALE targets. For all features XSX_S and h∈f^,fh∈\ f,f\, the theoretical uncentered ALE ALE~h(xS) ALE_h(x_S) (Eq. (2)) exists for all xSx_S. Moreover, the centered theoretical ALE ALEh(xS)ALE_h(x_S) exists and has finite second moment w.r.t. XSX_S, i.e., XS[ALEh(XS)2]<∞.E_X_S\! [ALE_h(X_S)^2 ]<∞. 0.A.2 Bias-Variance-Decomposition of the Estimators (Eq. (4)) Proof Fix a feature index S and an evaluation point xSx_S. Let D denote a random dataset used to estimate the feature effect. For better readability, we omit the point xSx_S. Under Assumption (i) for either dependent or independent D and f f, FD∣f^[(PDf−PD^f^)2]=FD∣f^[PDf2−2PDfPD^f^+PD^f^2] _FE_D f [(PD_f- PD_ f)^2 ]=E_FE_D f[PD_f^2-2PD_f PD_ f+ PD_ f^2] =PDf2−2PDfFD∣f^[PD^f^]+FD∣f^[PD^f^2] =PD_f^2-2PD_fE_FE_D f[ PD_ f]+E_FE_D f[ PD_ f^2] =PDf2−2PDfFD∣f^[PD^f^]+FVarD∣f^[PD^f^]+F[D∣f^[PD^f^]2] =PD_f^2-2PD_fE_FE_D f[ PD_ f]+E_FVar_D f[ PD_ f]+E_F[E_D f[ PD_ f]^2] =PDf2−2PDfFD∣f^[PD^f^]+FVarD∣f^[PD^f^]+VarFD∣f^[PD^f^] =PD_f^2-2PD_fE_FE_D f[ PD_ f]+E_FVar_D f[ PD_ f]+Var_FE_D f[ PD_ f] +F[D∣f^[PD^f^]]2 +E_F[E_D f[ PD_ f]]^2 =(PDf−FD∣f^[PD^f^])2+VarFD∣f^[PD^f^]+FVarD∣f^[PD^f^]. =(PD_f-E_FE_D f[ PD_ f])^2+Var_FE_D f[ PD_ f]+E_FVar_D f[ PD_ f]. The proof for ALE follows by replacing PDfPD_f with ALEfALE_f and PD^f PD_ f with ALE^f ALE_ f in the equations above. For this replacement argument, it suffices that all involved quantities are well-defined (Assumption (vi)) and have finite second moments (Assumption (v)), and that f, and hence ALEfALE_f, is non-random w.r.t. ℙFP_F. ∎ 0.A.3 Model Bias of the PD (Eq. (6)) Proof Fix xSx_S. Consider the product probability space (f^,XS¯)∼ℙF⊗ℙXS¯( f,X_ S) _F _X_ S (cf. Assumption(i)). By Assumption (i), f^(xS,XS¯) f(x_S,X_ S) is integrable under ℙF⊗ℙXS¯P_F _X_ S, hence PDf^(xS)=XS¯[f^(xS,XS¯)]PD_ f(x_S)=E_X_ S[ f(x_S,X_ S)] exists for ℙFP_F-almost all realizations of f f and is integrable w.r.t. ℙFP_F. Now F[PDf^(xS)]−PDf(xS)E_F\! [PD_ f(x_S) ]-PD_f(x_S) is equal to F[XS¯[f^(xS,XS¯)]]−XS¯[f(xS,XS¯)]=XS¯[F[f^(xS,XS¯)]−f(xS,XS¯)]. _F\! [E_X_ S\! [ f(x_S,X_ S) ] ]-E_X_ S\! [f(x_S,X_ S) ]=E_X_ S\! [E_F\! [ f(x_S,X_ S) ]-f(x_S,X_ S) ]. This follows by an application of Fubini’s theorem (exchange of integrals) since FXS¯[|f^(xS,XS¯)|]<∞E_FE_X_ S[| f(x_S,X_ S)|]<∞ by Assumption (i). ∎ 0.A.4 Model Variance of the PD (Eq. (7)) Proof As in §0.A.3, fix xSx_S and suppose Assumptions (i) and (i) hold. Then, the variance of the theoretical model PD exists and is: VarF[PDf^(xS)]=F[(XS¯[f^(xS,XS¯)]−FXS¯[f^(xS,XS¯)])2].Var_F [PD_ f(x_S) ]=E_F [ (E_X_ S[ f(x_S,X_ S)]-E_FE_X_ S[ f(x_S,X_ S)] )^2 ]. Applying Fubini’s theorem as FXS¯[|f^(xS,XS¯)|]<∞E_FE_X_ S[| f(x_S,X_ S)|]<∞ by Assumption (i), we have FXS¯[f^(xS,XS¯)]=XS¯F[f^(xS,XS¯)].E_FE_X_ S[ f(x_S,X_ S)]=E_X_ SE_F[ f(x_S,X_ S)]. By linearity of expectations, VarF[PDf^(xS)]=F[(XS¯[f^(xS,XS¯)−F[f^(xS,XS¯)]])2].Var_F [PD_ f(x_S) ]=E_F [ (E_X_ S [ f(x_S,X_ S)-E_F[ f(x_S,X_ S)] ] )^2 ]. Defining U:=f^(xS,XS¯)−F[f^(xS,XS¯)]U:= f(x_S,X_ S)-E_F[ f(x_S,X_ S)] and applying Jensen’s inequality to ϕ(t)=t2φ(t)=t^2 yields: (XS¯[U])2≤XS¯[U2]. (E_X_ S[U] )^2 _X_ S[U^2]. Therefore, VarF[PDf^(xS)]≤FXS¯[U2].Var_F [PD_ f(x_S) ] _FE_X_ S[U^2]. Since U2≥0U^2≥ 0, Tonelli’s theorem completes the proof: VarF[PDf^(xS)]≤FXS¯[U2]=XS¯F[U2]=XS¯VarF[f^(xS,XS¯)].Var_F [PD_ f(x_S) ] _FE_X_ S[U^2]=E_X_ SE_F[U^2]=E_X_ SVar_F [ f(x_S,X_ S) ]. ∎ 0.A.5 Estimation Variance of the PD (Eq. (8)) Proof Fix xSx_S and a model f f. Suppose Assumption (i) holds, and consider i.i.d. draws D=XS¯(i)i=1nD=\X_ S^(i)\_i=1^n from ℙXS¯∣f^P_X_ S f, conditional on f f. Define Y^S(i):=f^(xS,XS¯(i)) Y_S^(i):= f(x_S,X_ S^(i)), then by construction the Y^S(i) Y_S^(i) are i.i.d. conditional on f f, and by Assumption (i), we have XS¯(i)∣f^[(Y^S(i))2]<∞E_X_ S^(i) f\! [( Y_S^(i))^2 ]<∞ for ℙFP_F-almost all f f. Thereby VarD∣f^[PD^f^(xS)]=VarD∣f^[1n∑i=1nY^S(i)]=1n2VarD∣f^[∑i=1nY^S(i)] _D f [ PD_ f(x_S) ]=Var_D f\! [ 1n _i=1^n Y_S^(i) ]= 1n^2Var_D f\! [ _i=1^n Y_S^(i) ] =1n2(∑i=1nVarD∣f^[Y^S(i)])=1n2(nVarD∣f^[Y^S(1)])=1nVarD∣f^[f^(xS,XS¯)], = 1n^2 ( _i=1^nVar_D f[ Y_S^(i)] )= 1n^2 (n\,Var_D f[ Y_S^(1)] )= 1nVar_D f\! [ f(x_S,X_ S) ], with XS¯X_ S in the last step being drawn conditional on f f, so from D∣f D f. Taking expectation over the training randomness of f f gives the proof. ∎ 0.A.6 ALE Estimator Bias w.r.t. Binned Population ALE (§4.2) Proof Fix xSx_S, bins IS(k)k=1K\I_S(k)\_k=1^K, and a realization of f f. Suppose Assumptions (i) and (v) hold. Moreover, assume D=X(i)i=1nD=\X^(i)\_i=1^n (used for estimation) are i.i.d. draws from ℙXP_X, independent of f f, and assume nS(k)>0n_S(k)>0 for all bins k≤kS(xS)k≤ k_S(x_S). Let B=(B(1),…,B(n))B=(B^(1),…,B^(n)) be the bin assignments of the samples, where B(i)=kS(XS(i))B^(i)=k_S(X_S^(i)). Then, conditional on the full assignment vector B, i.e., on the product event ⋂i=1nXS(i)∈IS(B(i)) _i=1^n\X_S^(i)∈ I_S(B^(i))\, the samples remain independent; in particular, for each k the subcollection XS¯(i):B(i)=k\X_ S^(i):B^(i)=k\ is i.i.d. with XS¯(i)∼ℙXS¯∣XS∈IS(k)X_ S^(i) _X_ S X_S∈ I_S(k). Therefore, using the law of total expectation and D⟂f^D f, we obtain D∣f^[ALE~^f^(xS)]=BD∣B[∑k=1kS(xS)1nS(k)∑i:XS(i)∈IS(k)Δf^(k,XS¯(i))] _D f\! [ ALE_ f(x_S) ]=E_BE_D B\! [ _k=1^k_S(x_S) 1n_S(k) _i:X_S^(i)∈ I_S(k) _ f(k,X_ S^(i)) ] =B[∑k=1kS(xS)1nS(k)∑i:XS(i)∈IS(k)D∣B[Δf^(k,XS¯(i))]] =E_B\! [ _k=1^k_S(x_S) 1n_S(k) _i:X_S^(i)∈ I_S(k)E_D B\! [ _ f(k,X_ S^(i)) ] ] =B[∑k=1kS(xS)XS¯∣XS∈IS(k)[Δf^(k,XS¯)]]=B[ALE~f^K(xS)]=ALE~f^K(xS). =E_B\! [ _k=1^k_S(x_S)E_X_ S X_S∈ I_S(k)\! [ _ f(k,X_ S) ] ]=E_B\! [ ALE_ f^K(x_S) ]= ALE_ f^K(x_S). Thus, ALE~^f ALE_ f is unbiased for ALE~f^K ALE_ f^K: D∣f^[ALE~^f^(xS)]−ALE~f^K(xS)=0E_D f\! [ ALE_ f(x_S) ]- ALE_ f^K(x_S)=0. ∎ 0.A.7 Model Bias of the ALE (§4.2) Proof Fix xSx_S and let f^∼ℙF f _F. Suppose Assumptions (i), (i), (iv), and (v) hold, so (XS,XS¯)∼ℙX(X_S,X_ S) _X is an independent population draw, i.e., (XS,XS¯)⟂f^(X_S,X_ S) f, and ALE~f^(xS) ALE_ f(x_S) and ALE~f(xS) ALE_f(x_S) exist by Assumption (vi). For K∈ℕK , define gK(f^):=∑k=1kSK(xS)XS¯∣XS∈ISK(k)[Δf^,SK(k,XS¯)],g_K( f):= _k=1^k_S^K(x_S)E_X_ S X_S∈ I_S^K(k) [ _ f,S^K(k,X_ S) ], so that ALE~f^(xS)=limK→∞gK(f^) ALE_ f(x_S)= _K→∞g_K( f) by definition. Then, F[ALE~f^(xS)] _F\! [ ALE_ f(x_S) ] =∫ℱlimK→∞gK(f^)dℙF(f^)=limK→∞∫ℱgK(f^)dℙF(f^) = _F _K→∞g_K( f)\,dP_F( f)= _K→∞ _Fg_K( f)\,dP_F( f) =limK→∞∑k=1kSK(xS)∫ℱXS¯∣XS∈ISK(k)[Δf^,SK(k,XS¯)]dℙF(f^) = _K→∞ _k=1^k_S^K(x_S) _FE_X_ S X_S∈ I_S^K(k)\! [ _ f,S^K(k,X_ S) ]\,dP_F( f) =limK→∞∑k=1kSK(xS)XS¯∣XS∈ISK(k)[F[Δf^,SK(k,XS¯)]]. = _K→∞ _k=1^k_S^K(x_S)E_X_ S X_S∈ I_S^K(k)\! [E_F\! [ _ f,S^K(k,X_ S) ] ]. The second equality follows from the Dominated Convergence Theorem. Indeed, gK(f^)→ALE~f^(xS)g_K( f)→ ALE_ f(x_S) pointwise in f f by definition, and for every K, |gK(f^)| |g_K( f) | ≤∑k=1kSK(xS)XS¯∣XS∈ISK(k)[|Δf^,SK(k,XS¯)|]≤∑k=1kSK(xS)esssupS¯|Δf^,SK(k,S¯)| ≤ _k=1^k_S^K(x_S)E_X_ S X_S∈ I_S^K(k)\! [ | _ f,S^K(k,X_ S) | ]≤ _k=1^k_S^K(x_S) *ess\,sup_x_ S | _ f,S^K(k,x_ S) | ≤esssupS¯∑k=1kSK(xS)|f^(zk,SK,S¯)−f^(zk−1,SK,S¯)| ≤ *ess\,sup_x_ S\! _k=1^k_S^K(x_S)\!\! | f(z_k,S^K,x_ S)- f(z_k-1,S^K,x_ S) | ≤esssupS¯TV(t↦f^(t,S¯))≤V(f^), ≤ *ess\,sup_x_ STV\! (t f(t,x_ S) )≤ V( f), where V(f^)V( f) is the integrable bound from Assumption (iv). The last equality above follows from Fubini’s theorem, using Assumption (v) (integrability) for ℙXS¯∣XS∈ISK(k),f^=ℙXS¯∣XS∈ISK(k)P_X_ S X_S∈ I_S^K(k), f=P_X_ S X_S∈ I_S^K(k) since (XS,XS¯)⟂f^(X_S,X_ S) f. By linearity of expectation, F[ALE~f^(xS)]−ALE~f(xS)= _F\! [ ALE_ f(x_S) ]- ALE_f(x_S)= limK→∞∑k=1kSK(xS)XS¯∣XS∈ISK(k)[F[Δf^,SK(k,XS¯)]−Δf,SK(k,XS¯)]. _K→∞ _k=1^k_S^K(x_S)E_X_ S X_S∈ I_S^K(k)\! [E_F\! [ _ f,S^K(k,X_ S) ]- _f,S^K(k,X_ S) ]. Thus, this bias term vanishes if the finite differences of f f are unbiased w.r.t. those of f. One sufficient condition is that f f is a pointwise unbiased estimator of f. ∎ 0.A.8 Model Variance of the ALE (Eq. (11)) Proof Fix xSx_S, let f^∼ℙF f _F, and suppose Assumptions (i), (i), and (v) hold, so that as in §0.A.7, for each bin index k, all expectations w.r.t. ℙXS¯∣ℐkP_X_ S _k are independent of f f. For readability, we abbreviate the event XS∈IS(k)\X_S∈ I_S(k)\ by ℐkI_k. Then: VarF[ALE~f^K(xS)]=VarF(∑k=1kS(xS)XS¯∣ℐk[Δf^(k,XS¯)]) _F\! [ ALE_ f^K(x_S) ]=Var_F\! ( _k=1^k_S(x_S)E_X_ S _k [ _ f(k,X_ S) ] ) =F[(∑k=1kS(xS)XS¯∣ℐk[Δf^(k,XS¯)]−F[∑k=1kS(xS)XS¯∣ℐk[Δf^(k,XS¯)]])2] =E_F\! [ ( _k=1^k_S(x_S)E_X_ S _k [ _ f(k,X_ S) ]-E_F\! [ _k=1^k_S(x_S)E_X_ S _k [ _ f(k,X_ S) ] ] )^2 ] =F[(∑k=1kS(xS)XS¯∣ℐk[Δf^(k,XS¯)−F[Δf^(k,XS¯)]])2]. =E_F\! [ ( _k=1^k_S(x_S)E_X_ S _k [ _ f(k,X_ S)-E_F [ _ f(k,X_ S) ] ] )^2 ]. The last line follows from the linearity of expectations and exchanging FE_F and XS¯∣ℐkE_X_ S _k by an application of Fubini’s theorem. Indeed, using Assumption (v) for f^⟂(XS,XS¯) f (X_S,X_ S) ensures that the resulting joint integrand is integrable. Define Uk:=XS¯∣ℐk[Δf^(k,XS¯)−F[Δf^(k,XS¯)]],U_k:=E_X_ S _k\! [ _ f(k,X_ S)-E_F [ _ f(k,X_ S) ] ], and let =(U1,…,UkS(xS))⊤.U= (U_1,…,U_k_S(x_S) ) . By the Cauchy-Schwarz inequality, (∑k=1kS(xS)Uk)2=⟨,⟩2≤∥2∥2=kS(xS)∑k=1kS(xS)Uk2. ( _k=1^k_S(x_S)U_k )^2= 1,U ^2≤ 1 ^2\, ^2=k_S(x_S) _k=1^k_S(x_S)U_k^2. Moreover, F[kS(xS)∑k=1kS(xS)Uk2]=kS(xS)∑k=1kS(xS)F[Uk2]E_F\! [k_S(x_S)\! _k=1^k_S(x_S)U_k^2 ]=k_S(x_S) _k=1^k_S(x_S)E_F\! [U_k^2 ] by linearity of expectation. By Assumption (v), applying Jensen’s inequality to φ(t)=t2 (t)=t^2, we obtain Uk2≤XS¯∣ℐk[(Δf^(k,XS¯)−F[Δf^(k,XS¯)])2].U_k^2 _X_ S _k [ ( _ f(k,X_ S)-E_F[ _ f(k,X_ S)] )^2 ]. Putting everything together and applying Tonelli’s theorem (the integrand Uk2U_k^2 is nonnegative) gives the proof: VarF[ALE~f^K(xS)] _F\! [ ALE_ f^K(x_S) ] ≤kS(xS)∑k=1kS(xS)FXS¯∣ℐk[(Δf^(k,XS¯)−F[Δf^(k,XS¯)])2] ≤ k_S(x_S) _k=1^k_S(x_S)E_FE_X_ S _k\! [ ( _ f(k,X_ S)-E_F [ _ f(k,X_ S) ] )^2 ] =kS(xS)∑k=1kS(xS)XS¯∣ℐkF[(Δf^(k,XS¯)−F[Δf^(k,XS¯)])2] =k_S(x_S) _k=1^k_S(x_S)E_X_ S _kE_F\! [ ( _ f(k,X_ S)-E_F [ _ f(k,X_ S) ] )^2 ] =kS(xS)∑k=1kS(xS)XS¯∣ℐkVarF[Δf^(k,XS¯)]. =k_S(x_S) _k=1^k_S(x_S)E_X_ S _kVar_F\! [ _ f(k,X_ S) ]. □ 0.A.9 Estimation Variance of the ALE (Eq. (12)) Proof Fix xSx_S and bins IS(k)I_S(k), suppose Assumptions (i), (v) and (vi) hold, i.p. with D=X(i)i=1nD=\X^(i)\_i=1^n being i.i.d. draws from ℙX∣f^P_X f, conditional on f f, and assume nS(k)>0∀kn_S(k)>0\ ∀ k. For any fixed k, the variance of the finite differences σk2(f^):=VarXS¯∣XS∈IS(k),f^[Δf^(k,XS¯)] _k^2( f):=Var_X_ S X_S∈ I_S(k), f[ _ f(k,X_ S)] is finite by Assumption (v): σk2(f^)<∞ _k^2( f)<∞. Now, for each k, conditional on a model f f and on the full bin assignment vector B of all samples, the variables Δf^(k,XS¯(i)):B(i)=k\ _ f(k,X_ S^(i)):B^(i)=k\ (all those falling into IS(k)I_S(k)) are i.i.d. with variance σk2(f^) _k^2( f). Since also nS(k)n_S(k) and the events XS(i)∈IS(k)\X_S^(i)∈ I_S(k)\ are fixed conditional on B and f f (similar to §0.A.6), this gives us VarD∣f^,B[μ^k(f^)]=σk2(f^)/nS(k)Var_D f,B[ μ_k( f)]= _k^2( f)/n_S(k), where we define the unaccumulated value of the k-th bin as μ^k(f^):=1nS(k)∑i:XS(i)∈IS(k)Δf^(k,XS¯(i)). μ_k( f):= 1n_S(k) _i:X_S^(i)∈ I_S(k) _ f (k,X_ S^(i) ). Moreover, conditional on B and f f, different bins use disjoint subsets of the i.i.d. sample, hence CovD∣f^,B[μ^k(f^),μ^ℓ(f^)]=0Cov_D f,B[ μ_k( f), μ_ ( f)]=0 for k≠ℓk≠ . Altogether we get: FB∣f^VarD∣f^,B[ALE~^f^(xS)]=FB∣f^VarD∣f^,B[∑k≤kS(xS)μ^k(f^)] _FE_B fVar_D f,B\! [ ALE_ f(x_S) ]=E_FE_B fVar_D f,B [ _k≤ k_S(x_S) μ_k( f) ] =∑k≤kS(xS)FB∣f^VarD∣f^,B[μ^k(f^)]=∑k≤kS(xS)FB∣f^[1nS(k)σk2(f^)]. = _k≤ k_S(x_S)E_FE_B fVar_D f,B[ μ_k( f)]= _k≤ k_S(x_S)E_FE_B f [ 1n_S(k) _k^2( f) ]. □ Appendix 0.B Simulation Details & Results 0.B.1 Models, Hyperparameters, and Model Performances The XGBoost implementation from xgboost [chen_xgboost_2016], and the GAM implementation from pyGAM [serven_pygam_2018] were used. For the GAM, we considered the number of basis functions (∈[5,50]∈[5,50]) and the penalization control parameter (∈[0.001,1000]∈[0.001,1000], log-uniform) as tuning parameters. For XGBoost, we used the parameter spaces from [probst_tunability_2019] (confined to trees as base learners). Hyperparameter configurations for the overfitting models (OF) were carefully hand-picked to achieve strong performance on the training data while performing relatively poorly on holdout data. The optimal hyperparameters (OT) were selected by tuning the models on a separate data sample of size n (12501250 or 1000010000) for 200 trials using a Tree-structured Parzen Estimator (TPE) [bergstra2011algorithms] with MSE on separate holdout data (1000010000 samples for reliable performance estimates) as minimization objective. For both OF and OT models, the final selected hyperparameter configurations can be found on Github (link). We evaluate model performance across all repetitions on both the training and test data, with 1000010000 test samples to obtain reliable performance estimates. As intended, overfitting models show better training performance than the optimally tuned models, but underperform on test data and exhibit higher variance in their generalization performance. A linear regression model is used as a baseline and is outperformed by all OT models. These performance metrics can be found on GitHub (link). 0.B.2 Further Results: Bias-Variance-Analysis Table 6: Results for PD on Friedman1 averaged over 100 grid points. Bold numbers indicate the minimum per metric-feature-model-size-combination. Feature x_1 x_3 x_5 Metric MSE Bias Var MSE Bias Var MSE Bias Var n=1250n=1250 GAM_OF train 0.1184 0.0720 0.1173 0.1285 0.0639 0.1289 0.1210 0.0631 0.1212 val 0.1968 0.0929 0.1949 0.2084 0.0907 0.2073 0.2302 0.0955 0.2289 CV 0.1062 0.0697 0.1049 0.1094 0.0650 0.1090 0.1057 0.0652 0.1051 GAM_OT train 0.0166 0.0190 0.0168 0.0146 0.0238 0.0145 0.0145 0.0166 0.0147 val 0.0264 0.0266 0.0266 0.0187 0.0236 0.0188 0.0199 0.0221 0.0201 CV 0.0165 0.0182 0.0167 0.0144 0.0236 0.0143 0.0143 0.0166 0.0146 XGB_OF train 0.1361 0.2297 0.0862 0.4853 0.6748 0.0310 0.1142 0.2767 0.0390 val 0.1915 0.2885 0.1120 0.6402 0.7744 0.0419 0.1718 0.3097 0.0786 CV 0.1196 0.2742 0.0460 0.5429 0.7258 0.0167 0.1055 0.2853 0.0249 XGB_OT train 0.0327 0.1081 0.0217 0.0393 0.1511 0.0170 0.0208 0.0854 0.0140 val 0.0416 0.1206 0.0280 0.0509 0.1752 0.0209 0.0245 0.0925 0.0165 CV 0.0334 0.1227 0.0190 0.0439 0.1718 0.0149 0.0209 0.0941 0.0125 n=10000n=10000 GAM_OF train 0.0094 0.0178 0.0094 0.0086 0.0173 0.0086 0.0096 0.0192 0.0096 val 0.0137 0.0212 0.0137 0.0117 0.0211 0.0117 0.0129 0.0235 0.0128 CV 0.0095 0.0178 0.0095 0.0087 0.0173 0.0086 0.0097 0.0189 0.0096 GAM_OT train 0.0005 0.0082 0.0005 0.0004 0.0024 0.0004 0.0004 0.0035 0.0004 val 0.0013 0.0076 0.0013 0.0005 0.0021 0.0005 0.0005 0.0049 0.0005 CV 0.0005 0.0082 0.0005 0.0004 0.0024 0.0004 0.0004 0.0036 0.0004 XGB_OF train 0.0232 0.1194 0.0092 0.0675 0.2523 0.0040 0.0302 0.1594 0.0049 val 0.0326 0.1370 0.0143 0.0924 0.2949 0.0057 0.0386 0.1776 0.0073 CV 0.0235 0.1314 0.0065 0.0860 0.2886 0.0028 0.0344 0.1753 0.0038 XGB_OT train 0.0060 0.0266 0.0055 0.0060 0.0374 0.0047 0.0044 0.0274 0.0038 val 0.0078 0.0318 0.0070 0.0073 0.0428 0.0057 0.0055 0.0327 0.0046 CV 0.0052 0.0293 0.0045 0.0056 0.0427 0.0039 0.0041 0.0310 0.0033 Table 7: Results for ALE on Friedman1 averaged over 100 grid points. Bold numbers indicate the minimum per metric-feature-model-size-combination. Feature x_1 x_3 x_5 Metric MSE Bias Var MSE Bias Var MSE Bias Var n=1250n=1250 GAM_OF train 0.1280 0.0809 0.1258 0.1290 0.0629 0.1295 0.1246 0.0642 0.1248 val 0.4060 0.1738 0.3892 0.4199 0.2139 0.3875 0.4212 0.1265 0.4196 CV 0.1745 0.1855 0.1451 0.1497 0.1479 0.1324 0.1212 0.0951 0.1162 GAM_OT train 0.0206 0.0336 0.0201 0.0146 0.0235 0.0145 0.0144 0.0165 0.0146 val 0.0795 0.1668 0.0535 0.0403 0.1219 0.0263 0.0316 0.1086 0.0205 CV 0.0400 0.1449 0.0196 0.0273 0.1159 0.0144 0.0257 0.1147 0.0130 XGB_OF train 3.6571 1.6329 1.0248 1.2235 0.8780 0.4683 2.6986 1.4763 0.5371 val 4.8262 0.5670 4.6600 1.5224 0.7511 0.9914 1.9319 0.4747 1.7654 CV 0.8781 0.4926 0.6573 0.7944 0.7932 0.1710 0.5237 0.4235 0.3562 XGB_OT train 0.0591 0.1629 0.0337 0.0344 0.0961 0.0260 0.0482 0.1701 0.0199 val 0.1411 0.2250 0.0937 0.1314 0.2818 0.0538 0.0770 0.1795 0.0463 CV 0.0878 0.2428 0.0299 0.0920 0.2690 0.0203 0.0563 0.1989 0.0173 n=10000n=10000 GAM_OF train 0.0107 0.0272 0.0103 0.0087 0.0171 0.0087 0.0096 0.0191 0.0096 val 0.0177 0.0312 0.0173 0.0127 0.0214 0.0127 0.0135 0.0237 0.0133 CV 0.0107 0.0267 0.0103 0.0088 0.0171 0.0088 0.0097 0.0189 0.0097 GAM_OT train 0.0016 0.0239 0.0011 0.0004 0.0023 0.0004 0.0004 0.0036 0.0004 val 0.0036 0.0233 0.0031 0.0005 0.0023 0.0005 0.0005 0.0048 0.0005 CV 0.0016 0.0239 0.0011 0.0004 0.0024 0.0004 0.0004 0.0036 0.0004 XGB_OF train 0.6796 0.8044 0.0336 0.4915 0.6802 0.0299 0.8473 0.9057 0.0279 val 0.2098 0.1518 0.1932 0.1711 0.3057 0.0803 0.1157 0.1788 0.0866 CV 0.0485 0.1308 0.0324 0.1091 0.3024 0.0183 0.0499 0.1821 0.0173 XGB_OT train 0.0070 0.0252 0.0066 0.0055 0.0209 0.0052 0.0045 0.0164 0.0044 val 0.0166 0.0392 0.0155 0.0112 0.0483 0.0091 0.0075 0.0291 0.0069 CV 0.0064 0.0276 0.0058 0.0062 0.0436 0.0044 0.0045 0.0307 0.0037 Table 8: Results for PD on Feynman I.29.16 averaged over 100 grid points. Bold numbers indicate the minimum per metric-feature-model-size-combination. Feature x_1 θ _1 d_1 Metric MSE Bias Var MSE Bias Var MSE Bias Var n=1250n=1250 GAM_OF train 0.2916 0.1357 0.2826 0.0660 0.0482 0.0659 0.0693 0.0497 0.0692 val 0.4493 0.1492 0.4417 0.0988 0.0617 0.0982 0.1011 0.0528 0.1017 CV 0.1577 0.0749 0.1574 0.0611 0.0458 0.0610 0.0650 0.0485 0.0648 GAM_OT train 0.0318 0.0358 0.0316 0.0199 0.0281 0.0197 0.0179 0.0194 0.0182 val 0.0430 0.0370 0.0431 0.0255 0.0310 0.0254 0.0215 0.0224 0.0217 CV 0.0299 0.0344 0.0297 0.0186 0.0273 0.0185 0.0170 0.0193 0.0172 XGB_OF train 0.0681 0.0623 0.0664 0.0032 0.0209 0.0029 0.0030 0.0105 0.0030 val 0.0878 0.0715 0.0855 0.0049 0.0245 0.0044 0.0039 0.0178 0.0037 CV 0.0429 0.0408 0.0426 0.0021 0.0217 0.0016 0.0014 0.0044 0.0014 XGB_OT train 0.0316 0.0746 0.0270 0.0095 0.0229 0.0093 0.0074 0.0138 0.0074 val 0.0404 0.0813 0.0349 0.0117 0.0227 0.0116 0.0083 0.0162 0.0083 CV 0.0297 0.0865 0.0230 0.0079 0.0221 0.0077 0.0063 0.0129 0.0064 n=10000n=10000 GAM_OF train 0.0146 0.0221 0.0146 0.0072 0.0214 0.0070 0.0071 0.0160 0.0071 val 0.0192 0.0253 0.0192 0.0090 0.0217 0.0088 0.0090 0.0184 0.0089 CV 0.0142 0.0219 0.0142 0.0070 0.0212 0.0068 0.0069 0.0156 0.0068 GAM_OT train 0.0061 0.0135 0.0061 0.0036 0.0179 0.0034 0.0033 0.0123 0.0033 val 0.0075 0.0166 0.0075 0.0040 0.0156 0.0039 0.0039 0.0134 0.0039 CV 0.0058 0.0134 0.0058 0.0034 0.0178 0.0032 0.0032 0.0121 0.0032 XGB_OF train 0.0103 0.0198 0.0103 0.0014 0.0293 0.0006 0.0007 0.0080 0.0006 val 0.0124 0.0207 0.0124 0.0017 0.0293 0.0008 0.0007 0.0078 0.0007 CV 0.0066 0.0156 0.0065 0.0012 0.0290 0.0004 0.0004 0.0069 0.0004 XGB_OT train 0.0057 0.0181 0.0056 0.0028 0.0203 0.0024 0.0023 0.0098 0.0023 val 0.0072 0.0204 0.0070 0.0033 0.0202 0.0030 0.0026 0.0106 0.0025 CV 0.0045 0.0185 0.0043 0.0024 0.0204 0.0020 0.0018 0.0093 0.0018 Table 9: Results for ALE on Feynman I.29.16 averaged over 100 grid points. Bold numbers indicate the minimum per metric-feature-model-size-combination. Feature x_1 θ _1 d_1 Metric MSE Bias Var MSE Bias Var MSE Bias Var n=1250n=1250 GAM_OF train 0.3122 0.1226 0.3075 0.0701 0.0591 0.0689 0.0738 0.0506 0.0737 val 0.7270 0.2335 0.6956 0.2011 0.0703 0.2029 0.1729 0.0797 0.1723 CV 0.2393 0.2227 0.1963 0.0712 0.0503 0.0711 0.0786 0.0556 0.0781 GAM_OT train 0.0337 0.0405 0.0332 0.0214 0.0400 0.0205 0.0179 0.0195 0.0181 val 0.0827 0.1494 0.0625 0.0322 0.0437 0.0314 0.0220 0.0230 0.0222 CV 0.0604 0.1697 0.0327 0.0185 0.0409 0.0174 0.0146 0.0183 0.0148 XGB_OF train 0.5272 0.3080 0.4473 0.1359 0.0539 0.1376 0.1118 0.0677 0.1109 val 1.8846 0.2021 1.9073 0.1266 0.0616 0.1270 0.0653 0.0428 0.0656 CV 0.3804 0.2820 0.3113 0.0381 0.0520 0.0366 0.0338 0.0502 0.0323 XGB_OT train 0.0362 0.0421 0.0356 0.0186 0.0379 0.0177 0.0129 0.0176 0.0130 val 0.1037 0.1760 0.0752 0.0318 0.0277 0.0321 0.0189 0.0240 0.0189 CV 0.0790 0.2141 0.0344 0.0120 0.0397 0.0108 0.0076 0.0146 0.0076 n=10000n=10000 GAM_OF train 0.0149 0.0257 0.0147 0.0079 0.0304 0.0072 0.0072 0.0162 0.0072 val 0.0219 0.0294 0.0218 0.0114 0.0308 0.0108 0.0096 0.0190 0.0096 CV 0.0148 0.0255 0.0147 0.0077 0.0303 0.0070 0.0070 0.0156 0.0070 GAM_OT train 0.0063 0.0191 0.0062 0.0041 0.0279 0.0035 0.0033 0.0123 0.0033 val 0.0081 0.0230 0.0078 0.0052 0.0268 0.0046 0.0040 0.0134 0.0039 CV 0.0061 0.0191 0.0059 0.0040 0.0278 0.0033 0.0032 0.0120 0.0032 XGB_OF train 0.0470 0.0386 0.0471 0.0199 0.0541 0.0175 0.0105 0.0331 0.0098 val 0.2569 0.0855 0.2582 0.0197 0.0568 0.0171 0.0129 0.0196 0.0129 CV 0.0358 0.0354 0.0357 0.0059 0.0459 0.0039 0.0037 0.0178 0.0035 XGB_OT train 0.0069 0.0154 0.0069 0.0045 0.0241 0.0041 0.0029 0.0116 0.0029 val 0.0154 0.0212 0.0154 0.0087 0.0293 0.0081 0.0043 0.0132 0.0043 CV 0.0058 0.0256 0.0054 0.0043 0.0293 0.0035 0.0022 0.0098 0.0022 Table 10: Decomposition results of total variance VarTotVar_Tot into model VarModVar_Mod and estimation variance VarEstVar_Est on SimpleNormalCorrelated averaged over 100 grid points. Bold numbers indicate estimation strategy with minimal variance. Feature x_1 x_2 x_3 Metric VarTotVar_Tot VarModVar_Mod VarEstVar_Est VarTotVar_Tot VarModVar_Mod VarEstVar_Est VarTotVar_Tot VarModVar_Mod VarEstVar_Est PD XGB_OF train 0.0817 0.0813 0.0004 0.2263 0.2259 0.0004 0.0010 0.0009 0.0001 val 0.0857 0.0837 0.0020 0.2416 0.2395 0.0021 0.0014 0.0011 0.0003 CV 0.0505 0.0501 0.0004 0.1291 0.1287 0.0004 0.0007 0.0006 0.0001 XGB_OT train 0.0113 0.0112 0.0001 0.0141 0.0141 0.0001 0.0014 0.0014 0.0000 val 0.0123 0.0119 0.0004 0.0155 0.0151 0.0004 0.0019 0.0019 0.0000 CV 0.0090 0.0089 0.0001 0.0104 0.0103 0.0001 0.0013 0.0013 0.0000 ALE XGB_OF train 0.1012 0.0161 0.0851 0.1011 0.0249 0.0762 0.0741 0.0636 0.0105 val 0.5281 0.0359 0.4923 0.6017 0.1743 0.4275 0.0864 0.0424 0.0439 CV 0.1520 0.0620 0.0900 0.0994 0.0110 0.0884 0.0121 0.0003 0.0118 XGB_OT train 0.0088 0.0077 0.0011 0.0119 0.0107 0.0011 0.0066 0.0062 0.0005 val 0.0263 0.0124 0.0140 0.0275 0.0089 0.0186 0.0033 0.0003 0.0029 CV 0.0077 0.0049 0.0029 0.0103 0.0065 0.0038 0.0020 0.0014 0.0007 Table 11: Decomposition results of total variance VarTotVar_Tot into model VarModVar_Mod and estimation variance VarEstVar_Est on Feynman I.29.16 averaged over 100 grid points. Bold numbers indicate estimation strategy with minimal variance. Feature x_1 θ _1 d_1 Metric VarTotVar_Tot VarModVar_Mod VarEstVar_Est VarTotVar_Tot VarModVar_Mod VarEstVar_Est VarTotVar_Tot VarModVar_Mod VarEstVar_Est PD XGB_OF train 0.0664 0.0652 0.0013 0.0029 0.0027 0.0002 0.0030 0.0029 0.0002 val 0.0855 0.0794 0.0061 0.0044 0.0035 0.0009 0.0037 0.0029 0.0008 CV 0.0426 0.0414 0.0012 0.0016 0.0015 0.0002 0.0014 0.0013 0.0001 XGB_OT train 0.0270 0.0266 0.0004 0.0093 0.0091 0.0002 0.0074 0.0074 0.0000 val 0.0349 0.0334 0.0015 0.0116 0.0108 0.0007 0.0083 0.0081 0.0001 CV 0.0230 0.0226 0.0004 0.0077 0.0076 0.0002 0.0064 0.0064 0.0000 ALE XGB_OF train 0.4473 0.0498 0.3975 0.1376 0.1135 0.0240 0.1109 0.0864 0.0245 val 1.9073 -* 2.3566 0.1270 -* 0.1507 0.0656 -* 0.1164 CV 0.3113 -* 0.4167 0.0366 0.0099 0.0267 0.0323 0.0071 0.0252 XGB_OT train 0.0356 0.0302 0.0054 0.0177 0.0141 0.0037 0.0130 0.0112 0.0017 val 0.0752 0.0186 0.0565 0.0321 0.0096 0.0226 0.0189 0.0076 0.0113 CV 0.0344 0.0227 0.0116 0.0108 0.0065 0.0043 0.0076 0.0053 0.0024 * Due to instabilities in variance estimation for ALE, likely caused by unstable bin assignments, the estimated VarEstVar_Est sometimes exceeds the estimated VarTotVar_Tot. We omit VarModVar_Mod in these cases. (a) Mean estimation error PD (b) Mean estimation error ALE Figure 3: Mean estimation errors on SimpleNormalCorrelated for X1X_1 and X3X_3. For each sample size n, variances are averaged over all grid points. Axes are log-scale.