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.
Figures & tables
Figure 1: Conceptual overview of the four error components.
\begin{array}[]{l}X_{1},X_{2}\sim\mathrm{LogU}(0.1,10),\\ \theta_{1},\theta_{2}\sim\mathrm{U}(0,2\pi),\\ D_{1},D_{2}\sim\mathrm{U}(0,1),\text{ all indep.}\end{array}
Table 2: Data settings with underlying functions and feature structure.
Feature
x1
x2
x3
Metric
MSE
Bias
Var
MSE
Bias
Var
MSE
Bias
Var
n=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
Table 3: Results for PD on SimpleNormalCorrelated averaged over 100 grid points. Bold numbers are minimum per metric-feature-model-size-combination.
Feature
x1
x2
x3
Metric
MSE
Bias
Var
MSE
Bias
Var
MSE
Bias
Var
n=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
Table 4: Results for ALE on SimpleNormalCorrelated averaged over 100 grid points. Bold numbers are minimum per metric-feature-model-size-combination.
Feature
x1
x4
x7
Metric
VarTot
VarMod
VarEst
VarTot
VarMod
VarEst
VarTot
VarMod
VarEst
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
Table 5: Decomposition results of total variance VarTot into model variance VarMod and estimation variance VarEst on Friedman1 averaged over 100 grid points. Bold numbers indicate estimation strategy with minimal variance.
Figure 2: Mean estimation errors on Friedman1 for X1 and X4 . For each sample size n , the variances are averaged over all grid points. Both axes are log-scale.
Appendix figures & tables7 assets
Supplementary material from the paper’s appendix.
Appendix
Feature
x1
x3
x5
Metric
MSE
Bias
Var
MSE
Bias
Var
MSE
Bias
Var
n=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
Appendix
Table 6: Results for PD on Friedman1 averaged over 100 grid points. Bold numbers indicate the minimum per metric-feature-model-size-combination.
Feature
x1
x3
x5
Metric
MSE
Bias
Var
MSE
Bias
Var
MSE
Bias
Var
n=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
Appendix
Table 7: Results for ALE on Friedman1 averaged over 100 grid points. Bold numbers indicate the minimum per metric-feature-model-size-combination.
Feature
x1
θ1
d1
Metric
MSE
Bias
Var
MSE
Bias
Var
MSE
Bias
Var
n=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
Appendix
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
x1
θ1
d1
Metric
MSE
Bias
Var
MSE
Bias
Var
MSE
Bias
Var
n=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
Appendix
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
x1
x2
x3
Metric
VarTot
VarMod
VarEst
VarTot
VarMod
VarEst
VarTot
VarMod
VarEst
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
Appendix
Table 10: Decomposition results of total variance VarTot into model VarMod and estimation variance VarEst on SimpleNormalCorrelated averaged over 100 grid points. Bold numbers indicate estimation strategy with minimal variance.
Feature
x1
θ1
d1
Metric
VarTot
VarMod
VarEst
VarTot
VarMod
VarEst
VarTot
VarMod
VarEst
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
Appendix
Table 11: Decomposition results of total variance VarTot into model VarMod and estimation variance VarEst on Feynman I.29.16 averaged over 100 grid points. Bold numbers indicate estimation strategy with minimal variance.
Figure 3: Mean estimation errors on SimpleNormalCorrelated for X1 and X3 . For each sample size n , variances are averaged over all grid points. Axes are log-scale.
Statistical learning theory studies when data-driven algorithms generalize, how many samples they need, and what statistical or computational limits apply.
Feature importance methods using unrestricted permutations are flawed due to extrapolation errors; such errors appear in all non-trivial variable importance approaches. We propose three new approaches: conditional model reliance and Knockoffs with Gaussian transformation, and restricted ALE plot designs. Theoretical and numerical results show our strategies reduce/eliminate extrapolation.
Emanuele Borgonovo, Francesco Cappelli, Xuefei Lu +2
Department of Decision Sciences and Bocconi Institute for Data Science and Analytics, Bocconi University, Milan, Italy · Department of Decision Science, Bocconi University, Milan, Italy · SKEMA Business School, Université Côte d’Azur, France +2
Explanations of machine learning models are usually judged by criteria that are hard to compare. We propose a simpler test: if an explanation really describes how a model uses its features, it should be possible to rebuild the model's predictions from it. We turn each explanation into a predictor by reading each feature's effect and adding them up, and measure how well that predictor reproduces the model on unseen data. Nothing is fitted, so the score reflects the explanation itself. The test applies to any explanation that can be written as a function of the features; we demonstrate it on partial dependence plots (PDP), accumulated local effects (ALE), SHAP and LIME. We prove that summing partial dependence curves gives the best possible additive summary of a model when its features are independent, and that this fails when they are dependent. Across 13 real datasets and 9 synthetic designs and four model families, which method scores best depends entirely on feature dependence: where features are independent SHAP is slightly worse than PDP, exactly as the theory predicts; on dependent real data SHAP leads. Some widely used quality metrics even prefer a damaged explanation to an intact one.
Jacob Selbæk, Hugo L. Hammer
Department of Computer Science, Oslo Metropolitan University, Oslo, Norway · Department of Holistic Systems, SimulaMet, Oslo, Norway · Department of Plastic and Reconstructive Surgery, Oslo University Hospital, Oslo, Norway
We study fixed-design random feature ridge regression beyond the mean-kernel approximation. The expectation is taken over the frozen-feature ensemble, conditional on the training sample. Because the predictor is a nonlinear function of the empirical kernel, its mean training error, test error, and conditional generalization gap depend on centered kernel covariances as well as on mean kernel objects. We derive the covariance-level, or one-loop, corrections from a finite resolvent identity. This avoids an almost-sure Neumann-series assumption and yields an explicit remainder bound for the training error. The test correction additionally requires mixed train--test covariance tensors and therefore cannot, in general, be recovered from the train-restricted kernel law alone. Numerical experiments show that the covariance term removes the leading inverse-width discrepancy in a controlled regularization regime, while an ablation that discards the mixed train--test covariance fails for the test error even though the training correction is unchanged. They also identify a width--regularization boundary, measured by a resolvent-fluctuation diagnostic, beyond which the second-order truncation loses accuracy. All inverse-width statements are for fixed sample size, depth, design, and positive kernel-level regularization. The analysis concerns frozen features and does not model feature learning or neural-tangent-kernel evolution.
Taeyoung Kim
Center for AI and Natural Sciences, Korea Institute for Advanced Study, Seoul 02455, South Korea