diff --git a/include/maths/analytics/CTreeShapFeatureImportance.h b/include/maths/analytics/CTreeShapFeatureImportance.h index ce4952ab72..503c7c8551 100644 --- a/include/maths/analytics/CTreeShapFeatureImportance.h +++ b/include/maths/analytics/CTreeShapFeatureImportance.h @@ -73,7 +73,10 @@ class MATHS_ANALYTICS_EXPORT CTreeShapFeatureImportance { //! Compute inner node values as weighted average of the children (leaf) values. //! - //! The weights are the number of rows of \p frame reaching each node. + //! The weights are the number of rows of \p frame reaching each node. A node + //! below the root which no rows reach uses equal weights instead, matching the + //! even split used when computing SHAP values. If no rows reach the root its + //! value is NaN. static void computeInternalNodeValues(TTreeVec& forest); //! Get the maximum depth of any tree in \p forest. diff --git a/lib/maths/analytics/CBoostedTreeImpl.cc b/lib/maths/analytics/CBoostedTreeImpl.cc index 16c63d9344..78f5a28d8c 100644 --- a/lib/maths/analytics/CBoostedTreeImpl.cc +++ b/lib/maths/analytics/CBoostedTreeImpl.cc @@ -1517,9 +1517,16 @@ void CBoostedTreeImpl::initializeFixedCandidateSplits(core::CDataFrame& frame) { set.reserve(m_NumberSplitsPerFeature + 2); } for (auto row = beginRows; row != endRows; ++row) { + auto encodedRow = m_Encoder->encode(*row); for (std::size_t i = 0; i < features.size(); ++i) { if (state[i].size() <= m_NumberSplitsPerFeature + 1) { - state[i].insert(m_Encoder->encode(*row)[features[i]]); + // Missing values are NaN, which never compares equal to + // itself, so each one would otherwise count as a distinct + // value and corrupt the sorted candidate splits. + double value{encodedRow[features[i]]}; + if (core::CDataFrame::isMissing(value) == false) { + state[i].insert(value); + } } } } @@ -1538,7 +1545,10 @@ void CBoostedTreeImpl::initializeFixedCandidateSplits(core::CDataFrame& frame) { TDoubleVec values; for (std::size_t i = 0; i < features.size(); ++i) { - if (uniques[i].size() <= m_NumberSplitsPerFeature) { + // We need at least two distinct values to split. Since missing values are + // skipped there may be none, for example if the encoding was computed on a + // different data set. + if (uniques[i].size() > 1 && uniques[i].size() <= m_NumberSplitsPerFeature) { values.assign(uniques[i].begin(), uniques[i].end()); std::sort(values.begin(), values.end()); auto& featureCandidateSplits = m_FixedCandidateSplits[features[i]]; diff --git a/lib/maths/analytics/CTreeShapFeatureImportance.cc b/lib/maths/analytics/CTreeShapFeatureImportance.cc index 66eacea916..7dec3b28e1 100644 --- a/lib/maths/analytics/CTreeShapFeatureImportance.cc +++ b/lib/maths/analytics/CTreeShapFeatureImportance.cc @@ -169,6 +169,12 @@ void CTreeShapFeatureImportance::computeInternalNodeValues(TTree& tree, std::siz auto& rightChild = tree[node.rightChildIndex()]; double leftWeight{static_cast(leftChild.numberSamples())}; double rightWeight{static_cast(rightChild.numberSamples())}; + if (node.numberSamples() == 0 && nodeIndex != 0) { + // No counted rows reach this node so its parent gives it zero weight + // and its value doesn't change the expected value. Use the unweighted + // mean, which matches the even split shapRecursive uses for it. + leftWeight = rightWeight = 1.0; + } node.value((leftWeight * leftChild.value() + rightWeight * rightChild.value()) / (leftWeight + rightWeight)); } @@ -249,12 +255,21 @@ void CTreeShapFeatureImportance::shapRecursive(const TTree& tree, unwindPath(splitPath, pathIndex, nextIndex); } - double hotFractionZero{incomingFractionZero * - static_cast(tree[hotIndex].numberSamples()) / - static_cast(tree[nodeIndex].numberSamples())}; - double coldFractionZero{incomingFractionZero * - static_cast(tree[coldIndex].numberSamples()) / - static_cast(tree[nodeIndex].numberSamples())}; + // A node below the root which no counted rows reach leaves the child + // fractions as 0 / 0. An even split keeps them summing to one and doesn't + // change the expected value, since the fraction leading into the node is + // zero. It does affect the attributions of any row which some subset of its + // features routes into the node, for which it is the neutral choice. The + // root only has no samples if no rows were counted at all, which mustn't be + // masked. TreeInferenceModel in Elasticsearch must use the same rule. + double numberSamples{static_cast(tree[nodeIndex].numberSamples())}; + bool evenSplit{numberSamples == 0.0 && nodeIndex != 0}; + double hotFractionZero{ + incomingFractionZero * + (evenSplit ? 0.5 : static_cast(tree[hotIndex].numberSamples()) / numberSamples)}; + double coldFractionZero{ + incomingFractionZero * + (evenSplit ? 0.5 : static_cast(tree[coldIndex].numberSamples()) / numberSamples)}; this->shapRecursive(tree, encodedRow, hotIndex, hotFractionZero, incomingFractionOne, splitFeature, splitPath, nextIndex, shap); this->shapRecursive(tree, encodedRow, coldIndex, coldFractionZero, 0.0, diff --git a/lib/maths/analytics/unittest/CBoostedTreeTest.cc b/lib/maths/analytics/unittest/CBoostedTreeTest.cc index dbed0b6653..262e0a6540 100644 --- a/lib/maths/analytics/unittest/CBoostedTreeTest.cc +++ b/lib/maths/analytics/unittest/CBoostedTreeTest.cc @@ -62,6 +62,7 @@ class CBoostedTreeImplForTest { using TLossFunctionUPtr = maths::analytics::CBoostedTreeImpl::TLossFunctionUPtr; using TSizeVec = maths::analytics::CBoostedTreeImpl::TSizeVec; using TDoubleVec = maths::analytics::CBoostedTreeImpl::TDoubleVec; + using TFloatVecVec = maths::analytics::CBoostedTreeImpl::TFloatVecVec; using TBoostedTreeUPtr = std::unique_ptr; public: @@ -79,6 +80,10 @@ class CBoostedTreeImplForTest { return m_TreeImpl.featureSampleProbabilities(); } + const TFloatVecVec& fixedCandidateSplits() const { + return m_TreeImpl.m_FixedCandidateSplits; + } + void treeFeatureBag(TDoubleVec& probabilities, TSizeVec& treeFeatureBag) const { m_TreeImpl.treeFeatureBag(probabilities, treeFeatureBag); } @@ -914,6 +919,192 @@ BOOST_AUTO_TEST_CASE(testLowCardinalityFeatures) { BOOST_TEST_REQUIRE(rSquared > 0.94); } +BOOST_AUTO_TEST_CASE(testLowCardinalityFeaturesFixedCandidateSplitsWithMissingValues) { + + // Test that missing values don't leak into the fixed candidate splits of low + // cardinality features. Missing values are NaN, which is never equal to itself, + // so they must not be counted as extra distinct values. + + std::size_t rows{500}; + std::size_t cols{6}; + + test::CRandomNumbers rng; + TDoubleVecVec x(cols - 1); + for (std::size_t i = 0; i < cols - 1; ++i) { + rng.generateUniformSamples(0.0, 10.0, rows, x[i]); + for (auto& xj : x[i]) { + xj = std::floor(xj); + } + } + + // Make 2% of each feature's values missing. The target uses the complete values + // so that no training rows are dropped. + TDoubleVecVec xMissing{x}; + for (std::size_t i = 0; i < cols - 1; ++i) { + for (std::size_t j = i; j < rows; j += 50) { + xMissing[i][j] = core::CDataFrame::valueOfMissing(); + } + } + auto target = [&](const TRowRef& row) { + double result{0.0}; + for (std::size_t i = 0; i < cols - 1; ++i) { + result += static_cast(i + 1) * x[i][row.index()]; + } + return result; + }; + + auto frame = core::makeMainStorageDataFrame(cols, rows).first; + fillDataFrame(rows, 0, cols, xMissing, TDoubleVec(rows, 0.0), target, *frame); + + auto regression = maths::analytics::CBoostedTreeFactory::constructFromParameters( + 1, std::make_unique()) + .buildForTrain(*frame, cols - 1); + + // Each feature takes the values 0, 1, ..., 9 so its fixed candidate splits + // should be the midpoints 0.5, 1.5, ..., 8.5. + CBoostedTreeImplForTest treeImpl{regression->impl()}; + std::size_t numberFixed{0}; + for (const auto& splits : treeImpl.fixedCandidateSplits()) { + if (splits.empty() == false) { + ++numberFixed; + BOOST_REQUIRE_EQUAL(std::size_t{9}, splits.size()); + for (std::size_t j = 0; j < splits.size(); ++j) { + BOOST_REQUIRE_EQUAL(static_cast(j) + 0.5, + static_cast(splits[j])); + } + } + } + BOOST_REQUIRE_EQUAL(cols - 1, numberFixed); +} + +BOOST_AUTO_TEST_CASE(testLowCardinalityFeaturesWithMissingValues) { + + // Test training a linear model on low cardinality features with a few missing + // values. Corrupted candidate splits produce splits which no training rows reach + // and degrade the model. + + std::size_t trainRows{500}; + std::size_t testRows{200}; + std::size_t rows{trainRows + testRows}; + double noiseVariance{4.0}; + std::size_t cols{6}; + + test::CRandomNumbers rng; + TDoubleVecVec x(cols - 1); + for (std::size_t i = 0; i < cols - 1; ++i) { + rng.generateUniformSamples(0.0, 10.0, rows, x[i]); + for (auto& xj : x[i]) { + xj = std::floor(xj); + } + } + + // Make 2% of each feature's values missing. The target uses the complete values + // so that no training rows are dropped. + TDoubleVecVec xMissing{x}; + for (std::size_t i = 0; i < cols - 1; ++i) { + for (std::size_t j = i; j < rows; j += 50) { + xMissing[i][j] = core::CDataFrame::valueOfMissing(); + } + } + auto target = [&](const TRowRef& row) { + double result{0.0}; + for (std::size_t i = 0; i < cols - 1; ++i) { + result += static_cast(i + 1) * x[i][row.index()]; + } + return result; + }; + + TDoubleVec noise; + rng.generateNormalSamples(0.0, noiseVariance, rows, noise); + + auto frame = core::makeMainStorageDataFrame(cols, rows).first; + fillDataFrame(trainRows, testRows, cols, xMissing, noise, target, *frame); + + auto regression = maths::analytics::CBoostedTreeFactory::constructFromParameters( + 1, std::make_unique()) + .buildForTrain(*frame, cols - 1); + + regression->train(); + regression->predict(); + + // Every node should be reached by some of the rows its sample counts are + // computed from. + std::size_t numberZeroSampleNodes{0}; + for (const auto& tree : regression->trainedModel()) { + for (const auto& node : tree) { + if (node.numberSamples() == 0) { + ++numberZeroSampleNodes; + } + } + } + BOOST_REQUIRE_EQUAL(std::size_t{0}, numberZeroSampleNodes); + + double bias; + double rSquared; + std::tie(bias, rSquared) = computeEvaluationMetrics( + *frame, trainRows, rows, + [&](const TRowRef& row) { return regression->prediction(row)[0]; }, + target, noiseVariance / static_cast(rows)); + LOG_DEBUG(<< "bias = " << bias << ", rSquared = " << rSquared); + + BOOST_TEST_REQUIRE(rSquared > 0.9); +} + +BOOST_AUTO_TEST_CASE(testLowCardinalityFeatureAllMissingAfterSeparateEncoding) { + + // Test that we can train when a low cardinality feature selected by encoding + // on one data set has no values in the data set we train on. Skipping missing + // values then leaves no distinct values for its fixed candidate splits. + + std::size_t rows{300}; + std::size_t cols{4}; + + test::CRandomNumbers rng; + TDoubleVecVec x(cols - 1); + for (std::size_t i = 0; i < cols - 1; ++i) { + rng.generateUniformSamples(0.0, 10.0, rows, x[i]); + for (auto& xj : x[i]) { + xj = std::floor(xj); + } + } + auto target = [&](const TRowRef& row) { + double result{0.0}; + for (std::size_t i = 0; i < cols - 1; ++i) { + result += static_cast(i + 1) * x[i][row.index()]; + } + return result; + }; + + auto encodingFrame = core::makeMainStorageDataFrame(cols, rows).first; + fillDataFrame(rows, 0, cols, x, TDoubleVec(rows, 0.0), target, *encodingFrame); + + std::stringstream persistState; + { + auto encoded = maths::analytics::CBoostedTreeFactory::constructFromParameters( + 1, std::make_unique()) + .buildForEncode(*encodingFrame, cols - 1); + core::CJsonStatePersistInserter inserter(persistState); + encoded->acceptPersistInserter(inserter); + persistState.flush(); + } + + TDoubleVecVec xMissing{x}; + std::fill(xMissing[0].begin(), xMissing[0].end(), core::CDataFrame::valueOfMissing()); + auto trainingFrame = core::makeMainStorageDataFrame(cols, rows).first; + fillDataFrame(rows, 0, cols, xMissing, TDoubleVec(rows, 0.0), target, *trainingFrame); + + auto regression = maths::analytics::CBoostedTreeFactory::constructFromString(persistState) + .restoreFor(*trainingFrame, cols - 1); + regression->train(); + regression->predict(); + + trainingFrame->readRows(1, [&](const TRowItr& beginRows, const TRowItr& endRows) { + for (auto row = beginRows; row != endRows; ++row) { + BOOST_REQUIRE(std::isfinite(regression->prediction(*row)[0])); + } + }); +} + BOOST_AUTO_TEST_CASE(testLowTrainFractionPerFold) { // Test regression using a very low train fraction per fold. This should diff --git a/lib/maths/analytics/unittest/CTreeShapFeatureImportanceTest.cc b/lib/maths/analytics/unittest/CTreeShapFeatureImportanceTest.cc index 9499f1fde0..52fdad642b 100644 --- a/lib/maths/analytics/unittest/CTreeShapFeatureImportanceTest.cc +++ b/lib/maths/analytics/unittest/CTreeShapFeatureImportanceTest.cc @@ -24,6 +24,7 @@ #include #include +#include #include #include #include @@ -142,6 +143,59 @@ struct SFixtureSingleTree { mutable TTreeVec s_Trees; }; +struct SFixtureZeroSampleInnerNode { + SFixtureZeroSampleInnerNode() : s_Trees(1) { + + // Node 4 and its children have no samples. The counts are set directly so + // that rows can pass through node 4, as rows which weren't counted can when + // the model is used for inference. Each row of data takes a different path: + // it never reaches node 4, it reaches node 4's sibling, it passes through + // node 4 to either of its leaves, or it never reaches node 4 but its second + // feature alone routes it there. The last is the only case for rows which + // were counted. + TDoubleVecVec data{{0.7, 0.3}, {0.3, 0.3}, {0.3, 0.7}, {0.1, 0.7}, {0.7, 0.7}}; + + s_Frame = core::makeMainStorageDataFrame(s_NumberFeatures, s_NumberRows).first; + s_Frame->columnNames(columnNames(s_NumberFeatures)); + for (std::size_t i = 0; i < s_NumberRows; ++i) { + s_Frame->writeRow([&](core::CDataFrame::TFloatVecItr column, std::int32_t&) { + for (std::size_t j = 0; j < s_NumberFeatures; ++j, ++column) { + *column = data[i][j]; + } + }); + } + s_Frame->finishWritingRows(); + + CStubMakeDataFrameCategoryEncoder stubParameters{1, *s_Frame, 0}; + s_Encoder = std::make_unique(stubParameters); + + auto& tree = s_Trees[0]; + tree.resize(1); + tree[0].split(0, 0.5, true, 0.0, 0.0, 0.0, tree); + tree[1].split(1, 0.5, true, 0.0, 0.0, 0.0, tree); + tree[4].split(0, 0.2, true, 0.0, 0.0, 0.0, tree); + tree[2].value(toVector(20.0)); + tree[3].value(toVector(10.0)); + tree[5].value(toVector(7.0)); + tree[6].value(toVector(9.0)); + + TSizeVec numberSamples{4, 2, 2, 2, 0, 0, 0}; + for (std::size_t i = 0; i < tree.size(); ++i) { + tree[i].numberSamples(numberSamples[i]); + } + + s_TreeFeatureImportance = std::make_unique( + 1, *s_Frame, *s_Encoder, s_Trees, s_NumberFeatures); + } + + TDataFrameUPtr s_Frame; + std::size_t s_NumberFeatures{2}; + std::size_t s_NumberRows{5}; + TTreeShapFeatureImportanceUPtr s_TreeFeatureImportance; + TEncoderUPtr s_Encoder; + mutable TTreeVec s_Trees; +}; + struct SFixtureMultipleTrees { SFixtureMultipleTrees() : s_Trees(2) { @@ -481,6 +535,83 @@ BOOST_FIXTURE_TEST_CASE(testSingleTreeShap, SFixtureSingleTree) { }); } +BOOST_FIXTURE_TEST_CASE(testZeroSampleInnerNodeExpectedNodeValues, SFixtureZeroSampleInnerNode) { + + // A node with no samples has zero weight in its parent so it doesn't change the + // expected value. Its own value is the unweighted mean of its children, which + // matches the even split TreeSHAP uses for it. + TDoubleVec expectedValues{15.0, 10.0, 20.0, 10.0, 8.0, 7.0, 9.0}; + const auto& tree = s_Trees[0]; + for (std::size_t i = 0; i < tree.size(); ++i) { + BOOST_TEST_REQUIRE(tree[i].value()(0) == expectedValues[i]); + } + BOOST_TEST_REQUIRE(s_TreeFeatureImportance->baseline()(0) == 15.0); +} + +BOOST_FIXTURE_TEST_CASE(testZeroSampleInnerNodeShap, SFixtureZeroSampleInnerNode) { + + // The expected values are the Shapley values of the path-dependent value function + // v(S), computed by hand: v({}) = 15, v({f1}) = 20 or 10, v({f2}) = 15 for f2 = 0.3 + // and 0.5 * 20 + 0.5 * (0.5 * 7 + 0.5 * 9) = 14 for f2 = 0.7, and v({f1, f2}) is + // the prediction. For each row they sum to the prediction minus the baseline. + TDoubleVecVec expectedPhi{ + {5.0, 0.0}, {-5.0, 0.0}, {-5.0, -1.0}, {-6.0, -2.0}, {5.5, -0.5}}; + + s_Frame->readRows(1, [&](const TRowItr& beginRows, const TRowItr& endRows) { + for (auto row = beginRows; row != endRows; ++row) { + s_TreeFeatureImportance->shap( + *row, [&](const TSizeVec& indices, const TStrVec&, const TVectorVec& shap) { + BOOST_REQUIRE_EQUAL(indices.size(), row->numberColumns()); + for (auto i : indices) { + BOOST_REQUIRE_CLOSE_ABSOLUTE(expectedPhi[row->index()][i], + shap[i](0), 1e-7); + } + }); + } + }); +} + +BOOST_FIXTURE_TEST_CASE(testZeroSampleRootIsNotMasked, SFixtureZeroSampleInnerNode) { + + // The root only has no samples if no rows were counted at all. There is then no + // sample distribution to compute importance from, so this mustn't be masked by + // the even split used for zero-sample nodes below the root. + for (auto& node : s_Trees[0]) { + node.numberSamples(0); + } + maths::analytics::CTreeShapFeatureImportance treeFeatureImportance{ + 1, *s_Frame, *s_Encoder, s_Trees, s_NumberFeatures}; + + BOOST_TEST_REQUIRE(std::isnan(treeFeatureImportance.baseline()(0))); + s_Frame->readRows(1, [&](const TRowItr& beginRows, const TRowItr& endRows) { + for (auto row = beginRows; row != endRows; ++row) { + treeFeatureImportance.shap(*row, [&](const TSizeVec& indices, const TStrVec&, + const TVectorVec& shap) { + BOOST_REQUIRE_EQUAL(indices.size(), row->numberColumns()); + for (auto i : indices) { + BOOST_TEST_REQUIRE(std::isnan(shap[i](0))); + } + }); + } + }); +} + +BOOST_FIXTURE_TEST_CASE(testInconsistentSampleCountsAreNotMasked, SFixtureZeroSampleInnerNode) { + + // Counts computed from data always add up: a node's count is the sum of its + // children's. Here node 4 has samples but its children have none. Only a node + // which itself has no samples gets the even split, as in shapRecursive, so this + // mustn't produce a plausible looking baseline. + TSizeVec numberSamples{4, 2, 2, 0, 2, 0, 0}; + for (std::size_t i = 0; i < s_Trees[0].size(); ++i) { + s_Trees[0][i].numberSamples(numberSamples[i]); + } + maths::analytics::CTreeShapFeatureImportance treeFeatureImportance{ + 1, *s_Frame, *s_Encoder, s_Trees, s_NumberFeatures}; + + BOOST_TEST_REQUIRE(std::isnan(treeFeatureImportance.baseline()(0))); +} + BOOST_FIXTURE_TEST_CASE(testMultipleTreesShap, SFixtureMultipleTrees) { TStrVec expectedNames{s_Frame->columnNames()};