Source code for modeva.models.glmtree.glmtreeboost

import time

import numpy as np
from sklearn.base import BaseEstimator, RegressorMixin, ClassifierMixin
from sklearn.model_selection import train_test_split
from sklearn.utils import check_array
from sklearn.utils.extmath import softmax
from sklearn.utils.validation import check_is_fitted

from .glmtree import MoGLMTreeRegressor
from ..base import ModelBaseRegressor, ModelBaseClassifier
from ...auth import auth
from ...testsuite.interpret.fanova.base import InterpretFANOVA
from ...testsuite.interpret.fanova.tree_linear.extractor.glmtreeboost import InterpretBGLMTree
from ...utils.helper import limit_function_for_trial_license


class BaseGLMTreeBooster(BaseEstimator):

    def __init__(self,
                 n_estimators=100,
                 learning_rate=1.0,
                 n_epoch_no_change=5,
                 max_depth=1,
                 min_samples_leaf=50,
                 min_impurity_decrease=0,
                 split_custom=None,
                 n_screen_grid=1,
                 n_feature_search=5,
                 n_split_grid=20,
                 reg_lambda=None,
                 clip_predict=True,
                 verbose=True,
                 val_ratio=0.2,
                 random_state=0):

        auth.run()
        self.n_estimators = n_estimators
        self.learning_rate = learning_rate
        self.n_epoch_no_change = n_epoch_no_change

        self.max_depth = max_depth
        self.split_custom = split_custom
        self.min_samples_leaf = min_samples_leaf
        self.min_impurity_decrease = min_impurity_decrease
        self.n_screen_grid = n_screen_grid
        self.n_feature_search = n_feature_search
        self.n_split_grid = n_split_grid
        self.reg_lambda = reg_lambda
        self.clip_predict = clip_predict

        self.verbose = verbose
        self.val_ratio = val_ratio
        self.random_state = random_state

    def fit(self, X, y, sample_weight=None):
        """fit the GLMTree Boosting model

        Parameters
        ---------
        X : array-like of shape (n_samples, n_features)
            containing the input dataset
        y : array-like of shape (n_samples,)
            containing target values
        sample_weight : array-like of shape (n_samples,)
            containing the weight of each sample
        Returns
        -------
        object 
            self : Estimator instance.
        """

        start = time.time()
        X, y, sample_weight = self._validate_fit_inputs(X, y, sample_weight)

        self.estimators_ = []
        self.learning_rates_ = [1] + [self.learning_rate] * (self.n_estimators - 1)
        self.tr_idx_, self.val_idx_ = train_test_split(np.arange(X.shape[0]), test_size=self.val_ratio,
                                                       random_state=self.random_state)
        self._fit(X, y, sample_weight)
        self.time_cost_ = time.time() - start
        self.is_fitted_ = True
        return self

    def get_raw_output(self, x):
        """output

        Parameters
        ---------
        x : array-like of shape (n_samples, n_features)
            containing the input dataset
        Returns
        -------
        np.ndarray of shape (n_samples,),
        """
        pred = 0
        for indice, est in enumerate(self.estimators_):
            pred += self.learning_rates_[indice] * est.predict(x)
        return pred

    def extract_model_info(self, X, feature_names, feature_types):
        test_obj = InterpretBGLMTree(model=self,
                                     feature_names=feature_names,
                                     feature_types=feature_types,
                                     purification=True,
                                     is_boost=True)
        test_obj.fit(X)
        self.modeva_effects_ = test_obj.effects_
        self.modeva_intercept_ = test_obj.intercept_
        self.predict_effect = test_obj._predict_effect

    def interpret(self, dataset):

        if not dataset.is_built_in():
            limit_function_for_trial_license()
        self.extract_model_info(X=dataset.train_x,
                                feature_names=dataset.feature_names,
                                feature_types=dataset.feature_types)
        return InterpretFANOVA(model=self, dataset=dataset)


[docs] class MoGLMTreeBoostRegressor(RegressorMixin, BaseGLMTreeBooster, ModelBaseRegressor): """GLMTree Boosting regressor using residual-based boosting. A gradient boosting regressor that uses GLMTree base models to iteratively fit residuals. Each tree has both linear and non-linear components to capture complex relationships while maintaining interpretability. Parameters ---------- name : str, default=None Model identifier name for reference. n_estimators : int, default=100 Number of boosting rounds (trees) to fit. max_depth : int, default=1 Maximum tree depth. Model is most interpretable when depth=1. learning_rate : float, default=1.0 Shrinkage rate applied to each tree's contribution. n_epoch_no_change : int, default=5 Early stopping rounds - training stops if validation loss doesn't improve. min_samples_leaf : int, default=50 Minimum samples required in a leaf node. min_impurity_decrease : float, default=0 Minimum required decrease in impurity to split a node. split_custom : dict, default=None Custom split points specified per feature. n_screen_grid : int, default=1 Grid size for initial split point screening. n_feature_search : int, default=5 Number of features to consider after screening. n_split_grid : int, default=20 Grid size for fine-grained split point search. reg_lambda : float, default=0.1 L1 regularization strength parameter. clip_predict : bool, default=False Whether to clip predictions to training data range. simplified : bool, default=True Use simplified partial linear regression for splits. val_ratio : float, default=0.2 Proportion of data used for validation. verbose : bool, default=False Whether to print training progress. random_state : int, default=0 Random seed for reproducibility. Attributes ---------- estimators_ : list Fitted GLMTree models. n_features_in_ : int Number of input features. """ def __init__(self, name: str = None, n_estimators=100, max_depth=1, learning_rate=1.0, n_epoch_no_change=5, min_samples_leaf=50, min_impurity_decrease=0, split_custom=None, n_screen_grid=1, n_feature_search=5, n_split_grid=20, reg_lambda=0.1, clip_predict=False, simplified=True, verbose=False, val_ratio=0.2, random_state=0): self.name = name super().__init__(n_estimators=n_estimators, learning_rate=learning_rate, n_epoch_no_change=n_epoch_no_change, max_depth=max_depth, min_samples_leaf=min_samples_leaf, min_impurity_decrease=min_impurity_decrease, split_custom=split_custom, n_screen_grid=n_screen_grid, n_feature_search=n_feature_search, n_split_grid=n_split_grid, reg_lambda=reg_lambda, clip_predict=clip_predict, verbose=verbose, val_ratio=val_ratio, random_state=random_state) self.simplified = simplified def _more_tags(self): """ Internal function for skipping some sklearn estimator checks. """ return {"_xfail_checks": {"check_sample_weights_invariance": "zero sample_weight is not equivalent to removing samples"}} @staticmethod def _get_loss(label, pred, sample_weight=None): """method to calculate the MSE loss """ loss = np.average((label - pred) ** 2, axis=0, weights=sample_weight) return loss def _fit(self, X, y, sample_weight): """fit the GLMTree Boosting model Parameters --------- X : array-like of shape (n_samples, n_features) containing the input dataset y : array-like of shape (n_samples,) containing target values sample_weight : array-like of shape (n_samples,) containing the weight of each sample """ # Initialize the intercept z = y.copy() loss_opt = np.inf early_stop_count = 0 if self.verbose: print("#### MoGLMTreeBoost Training ####") for indice in range(self.n_estimators): estimator = MoGLMTreeRegressor(max_depth=self.max_depth, min_samples_leaf=self.min_samples_leaf, min_impurity_decrease=self.min_impurity_decrease, split_custom=self.split_custom, n_screen_grid=self.n_screen_grid, n_feature_search=self.n_feature_search, n_split_grid=self.n_split_grid, reg_lambda=self.reg_lambda, clip_predict=self.clip_predict, simplified=self.simplified, random_state=self.random_state) estimator.fit(X[self.tr_idx_], z[self.tr_idx_], sample_weight[self.tr_idx_]) z = z - self.learning_rates_[indice] * estimator.predict(X) # update loss_new = estimator._get_loss(z[self.val_idx_], np.zeros((len(self.val_idx_),)), sample_weight[self.val_idx_]) if self.verbose: print("Iteration %d with validation loss %0.5f" % (indice + 1, loss_new)) if loss_opt > loss_new: loss_opt = loss_new early_stop_count = 0 else: early_stop_count += 1 if early_stop_count > self.n_epoch_no_change: if self.verbose: print("Training is terminated as validation loss stops decreasing.") break self.estimators_.append(estimator) else: if self.verbose: print("""Training is terminated as max_epoch is reached.""") def _predict(self, X): """output prediction for given samples Parameters --------- X : array-like of shape (n_samples, n_features) containing the input dataset Returns ------- np.ndarray of shape (n_samples,) containing prediction """ X = check_array(X) check_is_fitted(self) pred = self.get_raw_output(X) return pred
[docs] class MoGLMTreeBoostClassifier(ClassifierMixin, BaseGLMTreeBooster, ModelBaseClassifier): """GLMTree Boosting classifier using residual-based boosting. A gradient boosting classifier that uses GLMTree base models to iteratively fit residuals. Each tree has both linear and non-linear components to capture complex relationships while maintaining interpretability. Parameters ---------- name : str, default=None Model identifier name for reference. n_estimators : int, default=100 Number of boosting rounds (trees) to fit. max_depth : int, default=1 Maximum tree depth. Model is most interpretable when depth=1. learning_rate : float, default=1.0 Shrinkage rate applied to each tree's contribution. n_epoch_no_change : int, default=5 Early stopping rounds - training stops if validation loss doesn't improve. min_samples_leaf : int, default=50 Minimum samples required in a leaf node. min_impurity_decrease : float, default=0 Minimum required decrease in impurity to split a node. split_custom : dict, default=None Custom split points specified per feature. n_screen_grid : int, default=1 Grid size for initial split point screening. n_feature_search : int, default=5 Number of features to consider after screening. n_split_grid : int, default=20 Grid size for fine-grained split point search. reg_lambda : float, default=0.1 L1 regularization strength parameter. clip_predict : bool, default=False Whether to clip predictions to training data range. simplified : bool, default=True Use simplified partial linear regression for splits. val_ratio : float, default=0.2 Proportion of data used for validation. verbose : bool, default=False Whether to print training progress. random_state : int, default=0 Random seed for reproducibility. Attributes ---------- estimators_ : list Fitted GLMTree models. n_features_in_ : int Number of input features. """ def __init__(self, name: str = None, n_estimators=100, learning_rate=1.0, n_epoch_no_change=5, max_depth=1, min_samples_leaf=50, min_impurity_decrease=0, split_custom=None, n_screen_grid=1, n_feature_search=5, n_split_grid=20, reg_lambda=0.1, clip_predict=False, simplified=True, verbose=False, val_ratio=0.2, random_state=0): self.name = name super().__init__(n_estimators=n_estimators, learning_rate=learning_rate, n_epoch_no_change=n_epoch_no_change, max_depth=max_depth, min_samples_leaf=min_samples_leaf, min_impurity_decrease=min_impurity_decrease, split_custom=split_custom, n_screen_grid=n_screen_grid, n_feature_search=n_feature_search, n_split_grid=n_split_grid, reg_lambda=reg_lambda, clip_predict=clip_predict, verbose=verbose, val_ratio=val_ratio, random_state=random_state) self.simplified = simplified def _more_tags(self): """ Internal function for skipping some sklearn estimator checks. """ return {"binary_only": True, "_xfail_checks": {"check_sample_weights_invariance": "zero sample_weight is not equivalent to removing samples"}} @staticmethod def _get_loss(label, pred, sample_weight=None): """method to calculate the cross entropy loss """ EPSILON = 1e-7 with np.errstate(divide="ignore", over="ignore"): pred = np.clip(pred, EPSILON, 1. - EPSILON) loss = - np.average(label * np.log(pred) + (1 - label) * np.log(1 - pred), axis=0, weights=sample_weight) return loss def _fit(self, X, y, sample_weight): """Fits the GLMTree boosting classifier model using gradient boosting. Implements gradient boosting training by iteratively fitting GLMTree base models to pseudo-residuals, using adaptive sample weights and early stopping based on validation loss. The model optimizes log loss for binary classification. Parameters ---------- X : array-like of shape (n_samples, n_features) Training data features matrix. y : array-like of shape (n_samples,) Binary target values (0 or 1). sample_weight : array-like of shape (n_samples,) Sample weights for weighted model fitting. If None, uniform weights are used. """ # Initialize the intercept z = y.copy() * 4 - 2 loss_opt = np.inf pred_train = 0 pred_val = 0 early_stop_count = 0 sw_adaptive = sample_weight.copy() for indice in range(self.n_estimators): estimator = MoGLMTreeRegressor(max_depth=self.max_depth, min_samples_leaf=self.min_samples_leaf, min_impurity_decrease=self.min_impurity_decrease, split_custom=self.split_custom, n_screen_grid=self.n_screen_grid, n_feature_search=self.n_feature_search, n_split_grid=self.n_split_grid, reg_lambda=self.reg_lambda, clip_predict=self.clip_predict, simplified=self.simplified, random_state=self.random_state) estimator.fit(X[self.tr_idx_], z[self.tr_idx_], sw_adaptive[self.tr_idx_]) pred_train += self.learning_rates_[indice] * estimator.predict(X[self.tr_idx_]) proba_train = 1 / (1 + np.exp(- np.clip(pred_train.ravel(), -8, 8))) pred_val += self.learning_rates_[indice] * estimator.predict(X[self.val_idx_]) proba_val = 1 / (1 + np.exp(- np.clip(pred_val.ravel(), -8, 8))) # update loss_new = self._get_loss(y[self.val_idx_].ravel(), proba_val, sample_weight[self.val_idx_]) if self.verbose: print("Iteration %d with validation loss %0.5f" % (indice + 1, loss_new)) if loss_opt > loss_new: loss_opt = loss_new early_stop_count = 0 else: early_stop_count += 1 if early_stop_count > self.n_epoch_no_change: if self.verbose: print("Early stop as validation loss does not decrease for certain iterations.") break self.estimators_.append(estimator) sw_adaptive[self.tr_idx_] = proba_train * (1 - proba_train) sw_adaptive[self.tr_idx_] /= np.sum(sw_adaptive[self.tr_idx_]) sw_adaptive[self.tr_idx_] = np.maximum(sw_adaptive[self.tr_idx_], 2 * np.finfo(np.float64).eps) with np.errstate(divide="ignore", over="ignore"): z[self.tr_idx_] = np.where(y.ravel()[self.tr_idx_], 1. / proba_train, -1. / (1. - proba_train)) z[self.val_idx_] = np.where(y.ravel()[self.val_idx_], 1. / proba_val, -1. / (1. - proba_val)) z = np.clip(z, a_min=-8, a_max=8) def _decision_function(self, X): """Computes raw decision scores for samples. Returns the raw decision values before sigmoid transformation, representing the model's confidence in predicting the positive class. Higher positive values indicate stronger confidence in class 1, while lower negative values indicate stronger confidence in class 0. Parameters ---------- X : array-like of shape (n_samples, n_features) Input samples for prediction. Returns ------- array-like of shape (n_samples,) Decision function values, where larger values indicate higher confidence in the positive class. """ X = check_array(X) check_is_fitted(self) pred = self.get_raw_output(X) return pred def _predict_proba(self, X): """Predicts class probabilities for samples. Computes probability estimates for each class by applying softmax transformation to the decision function outputs. Returns probabilities for both classes, where the probabilities sum to 1 for each sample. Parameters ---------- X : array-like of shape (n_samples, n_features) Input samples for prediction. Returns ------- array-like of shape (n_samples, 2) Probability estimates for each class, where [:, 0] contains probabilities for class 0 and [:, 1] contains probabilities for class 1. """ pred = self._decision_function(X) pred_proba = softmax(np.vstack([-pred, pred]).T / 2, copy=False) return pred_proba