Source-linked AI summary

Semiparametric Inference for Conditional Shapley Feature Importance

Agostino Gnasso

arXiv:2609.10313v3stat.MEstat.ML

TL;DR

Conditional Shapley importance addresses the lack of uncertainty quantification for dependence-aware feature attribution. The paper defines a global loss-based target and develops a one-step estimator with a squared-loss correction. Its reported results include nominal-coverage behavior and a test rejection rate of 0.05 under the null, alongside a bound on copula misspecification bias.

  • Problem

    Existing conditional Shapley methods largely provide point estimates, leaving valid large-sample inference for dependence-aware conditional importance unresolved.

    Source: §1 · p. 3
  • Method

    The paper defines a global conditional SAGE parameter and estimates it with a cross-fitted one-step procedure using a U-statistic-style squared-loss correction.

    Source: §1 · pp. 3–4
  • Results

    The Wald test’s rejection rate under the null is 0.05, and the squared-loss correction is reported to restore nominal coverage at moderate Monte Carlo batch size.

    Source: §1 · p. 4, §7 · p. 19
  • Takeaways & Limitations

    The conditional importance is a model-internal statement about feature use under the data law, not a causal query.

    Source: §7 · p. 18

    Wald intervals undercovered on correlated coordinates in the reported simulations.

    Source: §7 · p. 19

Abstract

from arXiv · show

Shapley values are widely used for post-hoc feature attribution, but most estimators return point quantities and do not quantify uncertainty, and popular implementations sample out-of-coalition features from their marginal distribution, which misattributes importance when features are dependent. This paper studies the conditional formulation, in which out-of-coalition features are integrated out under their true conditional distribution. The target is a global, loss-based importance that pairs a conditional value function with a SAGE-style loss aggregation. We propose a one-step estimator with K-fold cross-fitting and a U-statistic correction of the squared loss that removes the Monte Carlo bias of the naive plug-in; it is $\sqrt{n}$-consistent and asymptotically normal under double-robust rate conditions, and the resulting Wald interval attains nominal coverage. A Pinsker-type bound quantifies the bias from misspecifying the working copula class, while vine copulas keep conditional sampling tractable. In a Gaussian design study with n= 500, the empirical coverage of the 95% interval lies between 0.91 and 0.96 across all features, the test holds its Type-I rate at 0.05, and it reaches power one for moderate signals. Applied to the UCI Concrete and California Housing data, the method identifies the conditionally informative features with Bonferroni-controlled significance.

1 Introduction

This paper develops inference for conditional Shapley importance, addressing a gap in uncertainty quantification when features are dependent. It defines a global, loss-based parameter and an estimator with theoretical coverage guarantees and a bound on copula misspecification bias.

1 Introduction

  • Marginal sampling can evaluate predictions at out-of-distribution points and inflate the apparent importance of correlated variables.The conditional value function instead respects the feature dependence structure, though it requires estimating conditional distributions.Source: §1 · p. 2
  • The conditional formulation targets information about predictions under the data-generating law, whereas the marginal formulation addresses an interventional question; neither is uniformly preferable.The paper treats them as different estimands, with conditional attribution relevant to the former question.Source: §1 · p. 2
  • Prior conditional methods largely focused on point estimation, while existing inference addressed marginal targets; valid large-sample inference for conditional Shapley importance remained open.The conditional nuisance is a law over arbitrary coalitions, so existing rate conditions and orthogonality arguments do not transfer unchanged.Source: §1 · pp. 2–3
  • Conditional SAGE, Ψj(P), combines a conditional value function with loss-based aggregation to summarize feature j’s contribution to predictive variance.It is a global, non-negative scalar that marginalizes other features conditionally.Source: §1 · p. 3
  • The proposed one-step estimator uses K-fold cross-fitting and a U-statistic-style squared-loss estimate based on two independent Monte Carlo sub-batches.The correction removes the naive plug-in’s VarX−S|XS(f)/B bias and is reported to restore nominal coverage at moderate B; conditional-law draws require carrying the correction through the copula nuisance.Source: §1 · p. 4
  • The paper proves asymptotic normality and Wald coverage, and gives a Pinsker-type bound for bias from working-copula misspecification.The bound vanishes for nonparametric vines under the simplifying assumption.Source: §1 · p. 4

2 Background

The paper defines conditional Shapley attribution using the true conditional distribution of features outside each coalition, then aggregates it into a global loss-based importance. Copula representations support conditional sampling, and the paper extends influence-function inference to this conditional target.

2.1 Conditional Shapley values and the SAGE parameter

  • The conditional Shapley value integrates out features outside a coalition under their conditional law, rather than treating them independently of the coalition.The value function uses conditional means, with the empty coalition equal to E[f(X)] and the full coalition equal to f(x).Source: §2 · p. 5, §2.1 · p. 5
  • The conditional SAGE importance turns local conditional Shapley values into a global, loss-based summary of feature importance.Its loss-based value function is non-positive, equals minus Var(f(X)) for the empty coalition, and vanishes for the full feature set.Source: §2.1 · p. 5

2.2 Copulas, vines, and influence-function inference

  • The joint feature density decomposes into a copula density and marginal densities, providing a representation for conditional feature distributions.Vine copulas decompose the joint copula into bivariate conditional copulas, and sequential inversion of h-functions costs O(p^2) per draw.Source: §2.2 · p. 6
  • The paper extends one-step influence-function inference to the conditional Shapley functional, where the nuisance is the full joint copula.Prior variable-importance constructions used an independence-style sampling scheme.Source: §2.2 · p. 7

3 Methods

The methods correct the Monte Carlo bias in estimating squared conditional losses and use cross-fitting to construct a one-step estimator. The resulting sample-mean form supports Wald intervals without a separate plug-in correction term.

3.1 Influence function and debiased squared loss

  • The authors report that the influence-function simplification enables Wald intervals to attain nominal coverage at moderate sample sizes.The local-parameter alternative would require a density-weighted Riesz representer or kernel localization.Source: §3.1 · p. 7
  • The naive squared-loss plug-in is biased upward by the conditional variance divided by B, and the bias can be non-negligible at moderate B.Splitting conditional draws into independent sub-batches yields a degree-two U-statistic correction with conditional unbiasedness.Source: §3.1 · pp. 7–8

3.2 One-step cross-fitted estimator

  • Pointwise Shapley values are estimated using sampled feature permutations and conditional expectations under a fold-specific fitted copula.The stated implementation uses K = 3 folds and M = 30 permutations; its Gaussian working class is correctly specified for the Gaussian designs, while other settings probe misspecification.Source: §3.2 · pp. 8–9
  • The one-step estimator uses K-fold cross-fitting and has a sample-mean form, with a Wald interval formed from its estimated standard error.Under cross-fitting, the nuisance bias is second order and dominated by the central-limit term, so no separate plug-in correction term is needed.Source: §3.2 · p. 8

4 Theory

Under the stated regularity and nuisance-rate assumptions, the estimator is asymptotically normal and its Wald intervals attain nominal coverage. The theory also bounds the effect of copula misspecification and identifies the pseudo-true target in that setting.

4.1 Asymptotic normality and coverage

  • Under Assumptions 1–4 and K-fold cross-fitting, Theorem 1 establishes asymptotic normality; under Assumptions 1–4, Theorem 2 establishes nominal coverage, with convergence to N(0, 1) by Slutsky.Source: §4.1 · pp. 10–11
  • At the null boundary, the Wald limit is degenerate, although the test empirically holds size near 0.05 at n = 500.A score interval is available when boundary calibration is the primary concern.Source: §4.1 · p. 11
  • The one-step expansion also gives joint asymptotic normality across features, supporting simultaneous inference procedures.The paper notes Bonferroni control of family-wise error and a multiplier bootstrap for sharper simultaneous confidence regions.Source: §4.1 · p. 11

4.2 Robustness to copula misspecification

  • A Pinsker-type bound gives a quantitative check on the effect of copula misspecification, while intervals under a misspecified working class center on its pseudo-true parameter.That parameter is defined through the KL projection of the true copula onto the working class.Source: §4.2 · pp. 11–12

5 Simulations

The simulations assess interval calibration, testing power, sensitivity to dependence and copula misspecification, and the difference between conditional and marginal attribution. Results show strong calibration and power in baseline settings, with undercoverage near collinearity and a distinct conditional attribution for a dependent but unused feature.

Sim A: calibration

  • At n = 500, 95% Wald-interval coverage was 0.91–0.96 across all five features, with absolute bias below 0.013 on every feature.Signal-feature coverage was X1 0.95, X2 0.92, X3 0.92; null-feature coverage was X4 0.96, X5 0.91.Source: §5 · p. 13, Table 1 · p. 15

Sim B: power

  • At δ = 0, the rejection rate was 0.05; power reached 0.99 at δ = 0.25 and 1.00 for δ ≥ 0.5.The test uses nominal α = 0.05.Source: §5 · p. 13, Table 2 · p. 15

Sim D: copula misspecification

  • Under Gaussian-copula misspecification, coverage fell to 0.80 for X1 and 0.88 for X2, while independent features remained at ≥ 0.85.The localized bias accords with the bound: the copula KL deficit affects the conditional distribution of the misspecified pair.Source: §5 · pp. 13–14

Sim E: conditional versus marginal attribution

  • At n = 500 and R = 60, conditional importance for X2 is Ψ2 = 1.28 > 0 and flags it in all replications; marginal importance is zero and significant in only 5%.On the correlated pair, the conditional estimator undercovers (0.77 for X1 and 0.70 for X2).Source: §5 · p. 14

Sim C: coverage under correlation

  • At ρ = 0.9, coverage fell to 0.73 for X1 and 0.68 for X2, while independent features remained near nominal; underestimated standard errors explain the undercoverage.The point estimator remains consistent, but near-collinearity creates a finite-sample regularity problem.Source: §5 · p. 13, Table 3 · p. 15

6 Applications

The applications use the one-step SAGE-conditional estimator with a Gaussian-copula conditional sampler on Concrete and California Housing data. All eight features in each dataset survive Bonferroni correction, with the strongest reported contributions differing by application.

Application setup

  • The applications fit a gradient-boosting predictor and use Gaussian-copula conditional sampling with empirical marginals for the one-step SAGE-conditional estimator.The reported intervals are Wald 95% confidence intervals, with Bonferroni-adjusted p-values.Source: §6 · pp. 15–16

Concrete compressive strength

  • In Concrete (n = 800), all eight ingredients survive Bonferroni correction; age, cement, and water dominate at Ψ values of 108, 65, and 38 MPa2.Ash has a smaller contribution, Ψ = 5.25, consistent with its role as a partial cement substitute rather than an independent strength driver.Source: §6 · p. 16

California Housing

  • In California Housing (n = 2000), median income dominates by a factor of 2.9; all eight features survive Bonferroni at α = 0.05.Latitude and longitude jointly account for the geographic component, and the ordering agrees qualitatively with prior geo-economic analyses.Source: §6 · pp. 16–18

7 Discussion and Limitations

The discussion scopes the target as conditional attribution for a fixed or fitted prediction function, describing model-internal feature use rather than a causal query. It also identifies near-collinearity and high dimensionality as practical boundaries for inference and computation.

7 Discussion and Limitations

  • The target is conditional Shapley attribution for a fixed or fitted prediction function, a statement of model-internal feature use rather than a causal query.Causal interpretations require additional assumptions.Source: §7 · p. 18
  • For p ≳100, naive fitting is expensive because of O(p2) pair-copulas and an O(p3) per-evaluation factor.Truncated vines reduce fitting and sampling costs; for very high dimensions, the discussion recommends grouping correlated features.Source: §7 · pp. 18–19
  • The Wald test holds its nominal size at the null in Sim B, while score intervals are available for boundary calibration and Bonferroni, FDR, or multiplier-bootstrap procedures support simultaneous inference.The constraint Ψj ≥0 creates a boundary at the null; the bootstrap uses the joint Gaussian limit.Source: §7 · p. 19
  • Permutation sampling replaces the O(2p) exact Shapley sum, with Monte Carlo error O(M −1/2) that can be made negligible by choosing M ≳n.Conditional sampling via vine h-functions costs O(p2) per draw.Source: §7 · pp. 18–19
  • At |corr| ≳0.9, the Wald interval can undercover because near-singular conditioning inflates the nuisance-rate constant and makes the second-order remainder non-negligible.The parameter remains identified and the point estimator consistent, but the interval's undercoverage does not improve with n; grouped importance or a bootstrap variance is recommended.Source: §7 · pp. 19–20

Data Availability Statement

The article makes its datasets, reference implementation, and reproduction scripts publicly available.

Data Availability Statement

  • The Concrete Compressive Strength and California Housing datasets used in the article are publicly available.The sources are the UCI Machine Learning Repository and StatLib, respectively.Source: Data Availability Statement · p. 20
  • The shapCopula repository provides a reference implementation of the proposed estimator.It includes scripts that reproduce the article's simulation and application results.Source: Data Availability Statement · p. 20
  • The repository scripts reproduce the results reported in the article.The passage identifies the repository as openly available.Source: Data Availability Statement · p. 20

Supplementary Materials

The supplement provides proofs for Theorems 1 and 3 and additional numerical studies. The proof of Theorem 1 establishes asymptotic normality by separating a leading sample-mean term from a negligible nuisance remainder.

Supplementary Materials

  • The shapCopula software includes vine- and Gaussian-copula conditional samplers, Wald inference routines, and the article's public datasets.Reproducibility scripts cover each simulation design and the applications, reproducing the reported tables and figures.Source: Supplementary Materials · p. 21
  • The supplement contains proofs of Theorems 1 and 3, a tuning-parameter sensitivity analysis, an oracle-sampler ablation, and a near-collinearity study as sample size grows.Equation, theorem, assumption, and table numbers without an “S” prefix refer to the main text.Source: Supplementary Materials · p. 26

S1 Proof of Theorem 1

  • Under Assumptions 1–4, the proof decomposes the estimator error into a leading term and a remainder, then shows the scaled remainder is negligible.The nuisance rate condition is max_S ∥ˆv−k − vS∥L2 = oP(n−1/4).Source: S1 Proof of Theorem 1 · pp. 26–27
  • The leading term is a mean-zero i.i.d. sample mean with finite variance, so the classical central limit theorem yields asymptotic normality.Combining this limit with the negligible remainder gives √n(ˆΨj − Ψj) ⇒ N(0, σ²j).Source: S1 Proof of Theorem 1 · pp. 26–27

S2 Proof of Theorem 3

The proof compares the Shapley attribution under the true copula with that under its KL projection in the working class. It bounds the resulting difference using the shared marginals and total-variation control.

S2 Proof of Theorem 3

  • The proof compares the true copula with its KL projection onto the working class while keeping the feature marginals fixed.The Shapley value is linear in the value function, enabling a comparison through differences in that function.Source: S2 Proof of Theorem 3 · pp. 27–28
  • The value function is negative expected squared loss, and its discrepancy between copula models is bounded by comparing expectations under their feature laws.The argument uses boundedness of the prediction function and the value-function difference.Source: S2 Proof of Theorem 3 · p. 28
  • A Pinsker-type inequality bounds the total-variation terms using KL divergence, yielding a bound on attribution error from copula misspecification.The proof also invokes data processing to control the conditional distributions.Source: S2 Proof of Theorem 3 · p. 28

S3 Sensitivity to estimation choices

Across these supplementary checks, near-collinearity exposes a boundary for Wald coverage even as the point estimate converges. Oracle-sampler and tuning comparisons also show where finite-sample calibration is insensitive to copula fitting and estimator settings.

S3 Sensitivity to estimation choices

  • 95% Wald CI coverage stays in [0.86, 1.00] across the tested K, M, and B settings.Source: Table S1 · p. 29
  • K = 2→5 raises null-X4 coverage from 0.88 to 0.94 but increases runtime K-fold; M > 30 changes coverage little, while B = 20 is near-nominal and B = 40 adds no improvement.Source: S3 Sensitivity to estimation choices · p. 29

S4 Oracle-sampler ablation

  • Oracle and deployable Gaussian-copula results coincide within Monte Carlo error, so calibration is governed by one-step/cross-fitting rather than copula fit, as Theorem 1 predicts.The oracle sampler uses the exact Gaussian conditionals while leaving the rest of the estimator unchanged.Source: S4 Oracle-sampler ablation · p. 30, Table S2 · p. 30

S5 Behavior under near-collinearity as n grows

Loading 2609.10313v3…