References / explanations for contribution predictions

[david-cortes]

XGBoost offers an option pred_contribs, for which the docs explain to some degree what the result is supposed to contain, but they do not explain how those are computed, and do not cite any reference.

If I understand it correctly, those contributions are based on SHAP, which comes from this paper:

Lundberg, Scott M., and Su-In Lee.
"A unified approach to interpreting model predictions."
Advances in neural information processing systems 30 (2017).

Is that the case?

If I understand the paper correctly, when it computes the effect of a feature being excluded, it is supposed to take a weighted average of the predictions that would be obtained under the different possible values, but it appears XGBoost is using the sum of hessians, which would only correspond to a proportional split of the data for squared loss:

xgboost/src/predictor/interpretability/shap.cc

Line 57 in f416974

tree.Stat(tree.LeftChild(nidx)).sum_hess;

Where can I find some reference describing the methodology used here and these kinds of choices?

[trivialfis]

Is that the case?

Recently @RAMitchell has been working on https://arxiv.org/abs/2605.04497

[david-cortes]

@trivialfis Thanks. Could you add that reference to the relevant places in the documentation about contribution predictions? For example, in predict.xgb.Booster, predict.xgboost, and the online docs for parameters. Would also be helpful to have a tutorial article about it like there are for e.g. interaction constraints and similar.

@RAMitchell Curious about the detail with the hessians in the computation. The paper mentions:

Each edge e has a split feature f(e) and an edge weight we = nchild/nparent, the proportion of training samples flowing through that edge.

.. which is then used for:

Define the empty prediction for leaf v with leaf value val(v)

.. but the code uses the sum of hessians instead, regardless of the objective.

Is there some equivalency that I might be missing?

[RAMitchell]

Hessians in XGBoost are exactly equivalent to weight. From memory, I think this might be discussed in the original xgboost paper.

[david-cortes]

Hessians in XGBoost are exactly equivalent to weight.

They would be for objectives with a constant hessian, like reg:squarederror, but not for others like reg:logistic.

From memory, I think this might be discussed in the original xgboost paper.

From a look at it, the paper discusses a process similar to iteratively reweighted least squares, where hessians as used as weights in each subsequent approximation of the objective function inside a tree, but that wouldn't correspond with proportions of rows in the original data.

[david-cortes]

@RAMitchell Here's a concrete example where one can see a difference between the sum of hessians and the number of rows.

Here I'm using treelite to extract the info and to manage the nodes as Python objects.

Details

import xgboost as xgb
import treelite
from sklearn.datasets import fetch_california_housing
import numpy as np
import polars as pl
import json
import plotnine as p9

tmp = fetch_california_housing()
X = pl.DataFrame(tmp.data, schema=tmp.feature_names)
y = tmp.target

# Fits a single-tree model with differen objectives
model = xgb.train(
    dtrain=xgb.DMatrix(X, y),
    num_boost_round=1,
    params={
        "objective": "reg:squarederror",
        "seed": 123,
        "tree_method": "exact",
    }
)

model_gamma = xgb.train(
    dtrain=xgb.DMatrix(X, y),
    num_boost_round=1,
    params={
        "objective": "reg:gamma",
        "seed": 123,
        "tree_method": "exact",
    }
)

tl_model = treelite.frontend.from_xgboost(model)
tl_json = json.loads(tl_model.dump_as_json())

leaf_pred = model.predict(
    xgb.DMatrix(X),
    pred_leaf=True,
)
leaf_hess = [
    {"node_id": node["node_id"], "sum_hess": node["sum_hess"]}
    for node in tl_json["trees"][0]["nodes"]
]

tl_model_gamma = treelite.frontend.from_xgboost(model_gamma)
tl_json_gamma = json.loads(tl_model_gamma.dump_as_json())

leaf_pred_gamma = model_gamma.predict(
    xgb.DMatrix(X),
    pred_leaf=True,
)
leaf_hess_gamma = [
    {"node_id": node["node_id"], "sum_hess": node["sum_hess"]}
    for node in tl_json_gamma["trees"][0]["nodes"]
]

df_node_stats = (
    pl.DataFrame(leaf_hess)
    .with_columns(
        pl.col("node_id")
        .cast(pl.Int64)
    )
    .join(
        pl.DataFrame({"node_id": leaf_pred})
        .with_columns(
            pl.col("node_id")
            .cast(pl.Int64)
        )
        .group_by("node_id")
        .agg(pl.len().cast(pl.Float64).alias("num_rows"))
        ,
        on="node_id",
    )
)

df_node_stats_gamma = (
    pl.DataFrame(leaf_hess_gamma)
    .with_columns(
        pl.col("node_id")
        .cast(pl.Int64)
    )
    .join(
        pl.DataFrame({"node_id": leaf_pred_gamma})
        .with_columns(
            pl.col("node_id")
            .cast(pl.Int64)
        )
        .group_by("node_id")
        .agg(pl.len().cast(pl.Float64).alias("num_rows"))
        ,
        on="node_id",
    )
)

(
    p9.ggplot(
        pl.concat([
            df_node_stats
            .with_columns(obj=pl.lit("reg:squarederror"))
            ,
            df_node_stats_gamma
            .with_columns(obj=pl.lit("reg:gamma"))
        ])
        , p9.aes(x="num_rows", y="sum_hess")
    )
    + p9.theme_bw()
    + p9.geom_point(alpha=0.25)
    + p9.geom_smooth(alpha=0.5, color="navy", method="loess")
    + p9.facet_grid(cols="obj")
).show()

Result is as follows:

So as one can see from the plot, in the case of objectives with non-constant hessian, the sum of hessians do not necessarily reflect fractions of the data.

I also see that even for squared error, the contribution predictions do not match against exact shapley values if one were to obtain them by brute-forcing them, although they are quite close:

Details

def assign_nrows_to_json():
    curr_node = 0
    curr_X = tmp.data
    assign_nrows_to_json_recursive(curr_node, curr_X, tl_json["trees"][0]["nodes"])


def assign_nrows_to_json_recursive(
    curr_node: int,
    curr_X: np.ndarray,
    tl_json_nodes: dict,
):
    node_info = tl_json_nodes[curr_node]
    node_info["n_rows"] = curr_X.shape[0]
    if "leaf_value" in node_info:
        return
    
    val = curr_X[:, node_info["split_feature_id"]]
    thr = node_info["threshold"]
    if node_info["comparison_op"] == "<":
        mask_left = val < thr
    elif node_info["comparison_op"] == "<=":
        mask_left = val <= thr
    elif node_info["comparison_op"] == ">":
        mask_left = val > thr
    elif node_info["comparison_op"] == ">=":
        mask_left = val >= thr
    
    left_X = curr_X[mask_left]
    right_X = curr_X[~mask_left]
    assign_nrows_to_json_recursive(
        node_info["left_child"],
        left_X,
        tl_json_nodes,
    )
    assign_nrows_to_json_recursive(
        node_info["right_child"],
        right_X,
        tl_json_nodes,
    )

assign_nrows_to_json()

def predict_singletree(
    x: np.ndarray,
    tl_json: dict
) -> float:
    curr_node = 0
    return predict_singletree_recursive(
        x,
        tl_json["trees"][0]["nodes"],
        curr_node,
    ) + tl_json["base_scores"][0]

def predict_singletree_recursive(
    x: np.ndarray,
    tl_json_nodes: dict,
    curr_node: int
) -> float:
    node_info = tl_json_nodes[curr_node]
    if "leaf_value" in node_info:
        return node_info["leaf_value"]
    
    val = x[ node_info["split_feature_id"] ]
    if np.isnan(val):
        nrows_left = tl_json_nodes[ node_info["left_child"] ]["n_rows"]
        nrows_right = tl_json_nodes[ node_info["right_child"] ]["n_rows"]
        weight_left = nrows_left / (nrows_left + nrows_right)
        weight_right = nrows_right / (nrows_left + nrows_right)
        return (
            weight_left * predict_singletree_recursive(x, tl_json_nodes, node_info["left_child"])
            +
            weight_right * predict_singletree_recursive(x, tl_json_nodes, node_info["right_child"])
        )
    
    thr = node_info["threshold"]
    if node_info["comparison_op"] == "<":
        follow_left = val < thr
    elif node_info["comparison_op"] == "<=":
        follow_left = val <= thr
    elif node_info["comparison_op"] == ">":
        follow_left = val > thr
    elif node_info["comparison_op"] == ">=":
        follow_left = val >= thr
    if follow_left:
        return predict_singletree_recursive(x, tl_json_nodes, node_info["left_child"])
    else:
        return predict_singletree_recursive(x, tl_json_nodes, node_info["right_child"])

chosen_row = 20_000

np.testing.assert_allclose(
    predict_singletree(
        tmp.data[chosen_row],
        tl_json
    )
    ,
    model.predict(
        xgb.DMatrix(X[[chosen_row]]),
        output_margin=True,
    ),
    rtol=1e-6, atol=1e-6
)

def calc_brute_force_shapley(chosen_row):
    # Makes a grid of size 2^n with combinations of know and missing
    # value for each possible column
    ncols = tmp.data.shape[1]
    grid_vals = pl.DataFrame(
        {
            f"col{col}": [ tmp.data[chosen_row, col], None ]
            for col in range(ncols)
        }
    ).select(pl.all().implode())
    for col in range(ncols):
        grid_vals = grid_vals.explode(f"col{col}")
    grid_vals_np = grid_vals.to_numpy()

    pred_grid = [
        predict_singletree(
            grid_vals_np[row],
            tl_json,
        )
        for row in range(grid_vals_np.shape[0])
    ]
    grid_vals = (
        grid_vals
        .with_columns(
            pred=pl.Series(pred_grid)
        )
    )

    shapley_vals = np.empty(ncols)
    for col in range(ncols):
        shapley_vals[col] = (
            grid_vals
            .filter(
                pl.col(f"col{col}").is_not_null()
            )
            ["pred"].mean()
            -
            grid_vals
            .filter(
                pl.col(f"col{col}").is_null()
            )
            ["pred"].mean()
        )
    return shapley_vals

bruteforce_shapleys = calc_brute_force_shapley(chosen_row)

xgb_contribs = model.predict(
    xgb.DMatrix(X[[chosen_row]]),
    pred_contribs=True,
)[:, :-1]

print("bruteforce: ", bruteforce_shapleys)
print("xgb: ", xgb_contribs)

bruteforce:  [-0.2828531   0.00458511 -0.05503178  0.00029512  0.00038373 -0.03721721
 -0.02451519 -0.01924271]
xgb:  [[-0.28556395  0.00525349 -0.05760277  0.00035022  0.0003655  -0.03506732
  -0.02871053 -0.02379541]]

Am I missing something here?

BTW, in the case of gamma objective, the numbers would differ more:

bruteforce:  [-3.20209854e-01  8.54544766e-04 -9.96606701e-02  2.75180105e-05
  6.17688510e-05 -4.72511472e-02 -5.02451348e-02 -1.09154148e-01]
xgb:  [[-3.8357857e-01  1.2909863e-03 -1.2454319e-01  5.7190686e-05
   5.6303663e-05 -6.0839515e-02 -7.7852063e-02 -1.3111255e-01]]

[RAMitchell]

Boosting in effect reweights training instances (adaboost does this explicitly). So the probabilities here are with respect to the re-weightings of training data. If this is an advantage or disadvantage I don't know - but its convenient. You could run a background dataset through the tree and get counts in theory. I've seen arguments that this path dependant shap algorithm is true to the model and the interventional method is true to the data.

[david-cortes]

@RAMitchell A couple comments:

• Could you update the reference paper linked earlier (https://arxiv.org/abs/2605.04497) to document what XGBoost is actually computing? The paper currently mentions using number of observations for these calculations.
• I understand from the references cited in the paper that the result should match with shapley values, but in the example above, they differ by more than what can be expected from floating point error, even for squared loss which has constant hessians. Is the procedure in xgboost meant to calculate something else?
• It doesn't really make sense to me that the probabilities would be with respect to hessians instead of data frequencies, because then it would be theoretically possible for two identical trees to produce different contributions if they have the same splits and leaf values but different hessians, which would run counter to the intuitive idea of what the shapley values should calculate. Are you aware of some reference introducing this theoretical change of using sum of hessians instead of data fractions? Also @mayer79 any comment about using weights proportional to hessians instead of data frequencies when calculating expectations over omitted features in contributions?

[RAMitchell]

Yes I will update the paper when I get the chance.

I am not sure how you are getting your brute force shap values. Maybe you can post that code and I can look at it next week.

It doesn't really make sense to me that the probabilities would be with respect to hessians instead of data frequencies, because then it would be theoretically possible for two identical trees to produce different contributions if they have the same splits and leaf values but different hessians, which would run counter to the intuitive idea of what the shapley values should calculate.

This is a valid and key criticism of the path dependant shap algorithm. If I remember correctly, even if you do it with counts you can still make an inconsistent counter example.

We are going to overhaul the shap interfaces in xgboost and add some more options and docs but it might take a few months.

[david-cortes]

Maybe you can post that code and I can look at it next week.

Code is given above. It will show if you click on 'Details'. It follows the procedure as explained in the paper with proportional splitting according to data frequencies, evaluating all the 2^n combinations of subsets for a single observation, and shapley values calculated according to their theoretical definition.

[RAMitchell]

Ok I will check it when I get back but I think its a difference in your implementation somewhere.

There isn't really one theoretical shapley algorithm for trees. The shapley weights are fixed but differences come from how you interpret removing a feature from the set.

[david-cortes]

@mayer79 I see treeshap ends up using the data frequencies for SHAP calculations if passing argument recalculate=TRUE, but even there that parameter is set as default to FALSE for xgboost.unify. Any comments here about usage of hessians vs. frequencies?

[mayer79]

If I understand it correctly, those contributions are based on SHAP, which comes from this paper:

Lundberg, Scott M., and Su-In Lee.
"A unified approach to interpreting model predictions."
Advances in neural information processing systems 30 (2017).

Is that the case?

The 2017 paper is about model agnostic SHAP (Kernel SHAP, permutation SHAP etc). pred_contrib uses a different algo (TreeSHAP) with different values. TreeSHAP was also developed in 2017, released in 2018 (code) and finally published in the 2020 paper in Nature.

@article{lundberg2020,
  title={From local explanations to global understanding with explainable AI for trees},
  author={Lundberg, Scott M. and Erion, Gabriel and Chen, Hugh and DeGrave, Alex and Prutkin, Jordan M. and Nair, Bala and Katz, Ronit and Himmelfarb, Jonathan and Bansal, Nisha and Lee, Su-In},
  journal={Nature Machine Intelligence},
  volume={2},
  number={1},
  pages={2522-5839},
  year={2020},
  publisher={Nature Publishing Group}
}