Source code for modeva.models.glmtree.glmtree

from abc import abstractmethod
from warnings import simplefilter

import numpy as np
from sklearn.base import BaseEstimator, RegressorMixin, ClassifierMixin, is_regressor
from sklearn.exceptions import ConvergenceWarning
from sklearn.linear_model import LinearRegression, Lasso, LassoCV, LogisticRegression, LogisticRegressionCV
from sklearn.linear_model._base import LinearModel
from sklearn.linear_model._base import _preprocess_data
from sklearn.utils import check_array
from sklearn.utils.extmath import softmax
from sklearn.utils.validation import _deprecate_positional_args
from sklearn.utils.validation import check_is_fitted, _check_sample_weight

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

simplefilter("ignore", category=ConvergenceWarning)


class xLinearRegression(RegressorMixin, LinearModel):
    """
    Marginal Regression Linear Regression.

    Parameters
    ----------
    fit_intercept : bool, default=True
        Whether to calculate the intercept for this model. If set
        to False, no intercept will be used in calculations
        (i.e. data is expected to be centered).

    copy_X : bool, default=True
        If True, X will be copied; else, it may be overwritten.

    Attributes
    ----------
    coef_ : array of shape (n_features, ) or (n_targets, n_features)
        Estimated coefficients for the linear regression problem.
        If multiple targets are passed during the fit (y 2D), this
        is a 2D array of shape (n_targets, n_features), while if only
        one target is passed, this is a 1D array of length n_features.

    intercept_ : float or array of shape (n_targets,)
        Independent term in the linear model. Set to 0.0 if
        `fit_intercept = False`.
    """

    @_deprecate_positional_args
    def __init__(self, *, fit_intercept=True, copy_X=True):
        self.fit_intercept = fit_intercept
        self.copy_X = copy_X

    def fit(self, X, y, sample_weight=None):
        """
        Fit linear model.

        Parameters
        ----------
        X : {array-like, sparse matrix} of shape (n_samples, n_features)
            Training data

        y : array-like of shape (n_samples,)
            Target values. Will be cast to X's dtype if necessary

        sample_weight : array-like of shape (n_samples,), default=None
            Individual weights for each sample

        Returns
        -------
        self : returns an instance of self.
        """

        X, y = self._validate_data(X, y, y_numeric=True, multi_output=True)

        if sample_weight is not None:
            sample_weight = _check_sample_weight(sample_weight, X,
                                                 dtype=X.dtype)
            sample_weight = sample_weight / sample_weight.sum()

        X, y, X_offset, y_offset, X_scale = _preprocess_data(
            X, y, fit_intercept=self.fit_intercept,
            copy=self.copy_X, sample_weight=sample_weight)

        if sample_weight is not None:
            Xy = np.dot(X.T, y.ravel() * sample_weight)
            average = np.average(X, weights=sample_weight, axis=0)
            variance = np.average((X - average) ** 2, weights=sample_weight, axis=0) / (1 - (sample_weight ** 2).sum())
            self.coef_ = np.divide(Xy, variance, out=np.zeros_like(Xy), where=variance != 0)
            self.coef_ = self.coef_.T
        else:
            Xy = np.dot(X.T, y.reshape(-1, 1)) / X.shape[0]
            variance = np.var(X, axis=0, ddof=1).reshape(-1, 1)
            self.coef_ = np.divide(Xy, variance, out=np.zeros_like(Xy), where=variance != 0)
            self.coef_ = self.coef_.T

        if y.ndim == 1:
            self.coef_ = np.ravel(self.coef_)
        self._set_intercept(X_offset, y_offset, X_scale)
        return self


class GLMTree(BaseEstimator):
    """
        Base class for classification and regression.
     """

    def __init__(self,
                 base_estimator,
                 max_depth=3,
                 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=False,
                 simplified=True,
                 random_state=0):

        auth.run()
        self.base_estimator = base_estimator
        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.clip_predict = clip_predict
        self.reg_lambda = reg_lambda
        self.simplified = simplified
        self.random_state = random_state

    @abstractmethod
    def _build_root(self):
        pass

    @abstractmethod
    def _build_leaf(self, sample_indice):
        pass

    def _get_split_position(self, feature_indice, sorted_feature, sorted_indice, n_split_grid):

        EPSILON = 1e-7
        split_position = []
        if self.split_custom is not None:
            split_points = self.split_custom[feature_indice]
            for i, val in enumerate(split_points):
                loc = np.digitize(val, bins=sorted_feature, right=False)
                split_position.append(loc)
        else:
            feature_range = sorted_feature[-1] - sorted_feature[0]
            if feature_range >= EPSILON:
                idx = 0
                n_samples = len(sorted_indice)
                diff_flag = (sorted_feature[1:] - sorted_feature[:-1]) <= EPSILON
                for i, _ in enumerate(sorted_indice):

                    if i == (n_samples - 1):
                        continue

                    if ((i + 1) < self.min_samples_leaf) or ((n_samples - i - 1) < self.min_samples_leaf):
                        continue

                    if diff_flag[i]:
                        continue

                    percentage = (idx + 1) / (n_split_grid + 1)
                    if n_samples > (n_split_grid + 1) * self.min_samples_leaf:
                        if (i + 1) / n_samples < percentage:
                            continue
                    elif n_samples > 2 * self.min_samples_leaf:
                        if (i + 1 - self.min_samples_leaf) / (n_samples - 2 * self.min_samples_leaf) < percentage:
                            continue
                    elif (i + 1) != self.min_samples_leaf:
                        continue
                    idx += 1
                    split_position.append(i + 1)
            return split_position

    def _evaluate_splits_simplified(self, node_x, node_y, node_sw, sorted_indice, split_position):

        node_sw = node_sw.reshape(-1, 1)
        node_y = node_y.reshape(-1, 1)

        n_left = 0
        n_total = np.sum(node_sw)

        Xy_left = 0
        Xy_total = np.dot(node_x.T, node_y.ravel() * node_sw.ravel())

        sum_left = 0
        sum_total = np.sum(node_x * node_sw, 0)

        y_sum_left = 0
        y_sum_total = np.sum(node_y * node_sw)

        sq_sum_left = 0
        sq_node_x = node_x ** 2
        sq_sum_total = np.sum(sq_node_x * node_sw, 0)

        sq_sw_sum_left = 0
        sq_sw = node_sw ** 2
        sq_sw_sum_total = np.sum(sq_sw, 0)

        start_indice = 0
        left_impurity_list = []
        right_impurity_list = []
        for loc in split_position:
            end_indice = loc
            left_indice = sorted_indice[:end_indice]
            incremental_indice = sorted_indice[start_indice:end_indice]

            n_left += np.sum(node_sw[incremental_indice])
            sum_left += np.sum(node_x[incremental_indice] * node_sw[incremental_indice], 0)
            sq_sum_left += np.sum(sq_node_x[incremental_indice] * node_sw[incremental_indice], 0)
            y_sum_left += np.sum(node_y[incremental_indice] * node_sw[incremental_indice], 0)
            sq_sw_sum_left += np.sum(sq_sw[incremental_indice])
            Xy_left += np.sum(node_x[incremental_indice] * node_y[incremental_indice] * node_sw[incremental_indice], 0)
            X_offset_left = sum_left / n_left
            y_offset_left = y_sum_left / n_left

            correct_left = (n_left ** 2) / (n_left ** 2 - sq_sw_sum_left)
            var_left = (sq_sum_left / n_left - (sum_left / n_left) ** 2) * correct_left
            coef_left = np.divide(Xy_left / n_left - X_offset_left * y_offset_left, var_left,
                                  out=np.zeros_like(Xy_left), where=var_left != 0)
            intercept_left = y_offset_left - np.dot(X_offset_left, coef_left.T)

            left_estimator = LinearRegression()
            left_estimator.coef_ = coef_left
            left_estimator.intercept_ = intercept_left
            left_impurity = self._evaluate_estimator(left_estimator, node_x[left_indice],
                                                     node_y[left_indice].ravel(), node_sw[left_indice].ravel())

            right_indice = sorted_indice[end_indice:]

            n_right = n_total - n_left
            sum_right = sum_total - sum_left
            sq_sum_right = sq_sum_total - sq_sum_left
            y_sum_right = y_sum_total - y_sum_left
            sq_sw_sum_right = sq_sw_sum_total - sq_sw_sum_left
            Xy_right = Xy_total - Xy_left
            X_offset_right = sum_right / n_right
            y_offset_right = y_sum_right / n_right

            correct_right = (n_right ** 2) / (n_right ** 2 - sq_sw_sum_right)
            var_right = (sq_sum_right / n_right - (sum_right / n_right) ** 2) * correct_right
            coef_right = np.divide(Xy_right / n_right - X_offset_right * y_offset_right, var_right,
                                   out=np.zeros_like(Xy_right), where=var_right != 0)
            intercept_right = y_offset_right - np.dot(X_offset_right, coef_right.T)

            right_estimator = LinearRegression()
            right_estimator.coef_ = coef_right
            right_estimator.intercept_ = intercept_right
            right_impurity = self._evaluate_estimator(right_estimator, node_x[right_indice],
                                                      node_y[right_indice].ravel(), node_sw[right_indice].ravel())
            start_indice = end_indice

            left_impurity_list.append(left_impurity)
            right_impurity_list.append(right_impurity)
        return left_impurity_list, right_impurity_list

    def _evaluate_splits_estimator(self, node_x, node_y, node_sw, sorted_indice, split_position):

        EPSILON = 1e-7
        left_impurity_list = []
        right_impurity_list = []
        for loc in split_position:
            end_indice = loc
            left_indice = sorted_indice[:end_indice]
            try:
                self.base_estimator.fit(node_x[left_indice], node_y[left_indice], node_sw[left_indice])
            except:
                proba = np.clip(np.mean((node_y[left_indice],)), EPSILON, 1 - EPSILON)
                self.base_estimator.coef_ = np.zeros((1, self.n_features_in_))
                self.base_estimator.intercept_ = np.log(proba / (1 - proba))

            left_impurity = self._evaluate_estimator(self.base_estimator, node_x[left_indice],
                                                     node_y[left_indice].ravel(), node_sw[left_indice])

            right_indice = sorted_indice[end_indice:]
            try:
                self.base_estimator.fit(node_x[right_indice], node_y[right_indice], node_sw[right_indice])
            except:
                proba = np.clip(np.mean((node_y[left_indice],)), EPSILON, 1 - EPSILON)
                self.base_estimator.coef_ = np.zeros((1, self.n_features_in_))
                self.base_estimator.intercept_ = np.log(proba / (1 - proba))

            right_impurity = self._evaluate_estimator(self.base_estimator, node_x[right_indice],
                                                      node_y[right_indice].ravel(), node_sw[right_indice])

            left_impurity_list.append(left_impurity)
            right_impurity_list.append(right_impurity)
        return left_impurity_list, right_impurity_list

    def _evaluate_split(self, node_x, node_y, node_sw, n_split_grid, feature_indice):

        EPSILON = 1e-7
        best_feature = None
        best_position = None
        best_threshold = None
        best_impurity = np.inf
        best_left_impurity = np.inf
        best_right_impurity = np.inf

        current_feature = node_x[:, feature_indice]
        sorted_indice = np.argsort(current_feature)
        sorted_feature = current_feature[sorted_indice]
        split_position = self._get_split_position(feature_indice, sorted_feature, sorted_indice, n_split_grid)

        if self.simplified and is_regressor(self):
            left_impurity_list, right_impurity_list = self._evaluate_splits_simplified(node_x, node_y, node_sw,
                                                                                       sorted_indice, split_position)
        else:
            left_impurity_list, right_impurity_list = self._evaluate_splits_estimator(node_x, node_y, node_sw,
                                                                                      sorted_indice, split_position)

        for i, loc in enumerate(split_position):
            end_indice = loc

            left_impurity = left_impurity_list[i]
            right_impurity = right_impurity_list[i]
            left_indice = sorted_indice[:end_indice]
            right_indice = sorted_indice[end_indice:]
            current_impurity = (node_sw[left_indice].sum() * left_impurity + node_sw[
                right_indice].sum() * right_impurity) / node_sw.sum()

            if current_impurity < (best_impurity - EPSILON):
                best_position = loc
                best_feature = feature_indice
                best_impurity = current_impurity
                best_left_impurity = left_impurity
                best_right_impurity = right_impurity
                best_threshold = (sorted_feature[loc - 1] + sorted_feature[loc]) / 2

        res = np.hstack(
            [best_position, best_feature, best_impurity, best_left_impurity, best_right_impurity, best_threshold])
        return res

    def _screen_features(self, sample_indice):

        node_x = self.x_[sample_indice]
        node_y = self.y_[sample_indice]
        node_sw = self.sample_weight_[sample_indice]
        allres = [self._evaluate_split(node_x, node_y, node_sw, self.n_screen_grid, feature_indice)
                  for feature_indice in self.split_features_]
        stat = np.vstack(allres)
        feature_impurity = stat[:, 2]  # impurity
        split_feature_indices = np.argsort(feature_impurity)[:self.n_feature_search]
        important_split_features = np.array(self.split_features_)[split_feature_indices]
        return important_split_features

    def _node_split(self, sample_indice):

        node_x = self.x_[sample_indice]
        node_y = self.y_[sample_indice]
        node_sw = self.sample_weight_[sample_indice]
        allres = [self._evaluate_split(node_x, node_y, node_sw, self.n_split_grid, feature_indice)
                  for feature_indice in self.important_split_features_]

        stat = np.vstack(allres).astype("float64")
        if np.isnan(stat).any():
            node = {"feature": None, "threshold": None, "left": None, "right": None,
                    "impurity": np.inf, "left_impurity": np.inf, "right_impurity": np.inf}
        else:
            best_stat = stat[np.nanargmin(stat[:, 2])]
            best_position = int(best_stat[0])
            best_feature = int(best_stat[1])
            best_impurity = best_stat[2]
            best_left_impurity = best_stat[3]
            best_right_impurity = best_stat[4]
            best_threshold = best_stat[5]
            if best_position is not None:
                sorted_indice = np.argsort(node_x[:, best_feature])
                best_left_indice = sample_indice[sorted_indice[:best_position]]
                best_right_indice = sample_indice[sorted_indice[best_position:]]
            node = {"feature": best_feature, "threshold": best_threshold, "left": best_left_indice,
                    "right": best_right_indice,
                    "impurity": best_impurity, "left_impurity": best_left_impurity,
                    "right_impurity": best_right_impurity}
        return node

    def _add_node(self, parent_id, is_left, is_leaf, depth, feature, threshold, impurity, sample_indice):

        self.node_count_ += 1
        if parent_id is not None:
            if is_left:
                self.tree_[parent_id].update({"left_child_id": self.node_count_})
            else:
                self.tree_[parent_id].update({"right_child_id": self.node_count_})

        node_id = self.node_count_
        n_samples = len(sample_indice)
        if is_leaf:
            predict_func, estimator, best_impurity = self._build_leaf(sample_indice)
            node = {"node_id": node_id, "parent_id": parent_id, "depth": depth, "feature": feature,
                    "impurity": best_impurity,
                    "n_samples": n_samples, "is_left": is_left, "is_leaf": is_leaf,
                    "value": np.mean(self.y_[sample_indice]),
                    "predict_func": predict_func, "estimator": estimator}
            self.leaf_estimators_.update({node_id: estimator})
        else:
            node = {"node_id": node_id, "parent_id": parent_id, "depth": depth, "feature": feature,
                    "impurity": impurity,
                    "n_samples": n_samples, "is_left": is_left, "is_leaf": is_leaf,
                    "value": np.mean(self.y_[sample_indice]),
                    "left_child_id": None, "right_child_id": None, "threshold": threshold}
        self.tree_.update({node_id: node})
        return node_id

    def fit(self, X, y, sample_weight=None):

        # initialize
        EPSILON = 1e-7
        self.tree_ = {}
        self.node_count_ = 0
        self.leaf_estimators_ = {}
        self.x_, self.y_, self.sample_weight_ = self._validate_fit_inputs(X, y, sample_weight)
        sample_indice = np.arange(self.x_.shape[0])

        if self.reg_lambda is None:
            self.reg_lambda_ = [0.0]
        else:
            self.reg_lambda_ = self.reg_lambda
        if self.split_custom is None:
            self.split_features_ = np.arange(self.n_features_in_).tolist()
            if self.n_feature_search > self.n_features_in_ and (self.max_depth >= 1):
                self.important_split_features_ = self.split_features_
            else:
                self.important_split_features_ = self._screen_features(sample_indice)
        else:
            self.important_split_features_ = list(self.split_custom.keys())

        np.random.seed(self.random_state)
        root_impurity = self._build_root()
        root_node = {"sample_indice": sample_indice,
                     "parent_id": None,
                     "depth": 0,
                     "impurity": root_impurity,
                     "is_left": False}
        pending_node_list = [root_node]
        while len(pending_node_list) > 0:
            stack_record = pending_node_list.pop()
            sample_indice = stack_record["sample_indice"]
            parent_id = stack_record["parent_id"]
            depth = stack_record["depth"]
            impurity = stack_record["impurity"]
            is_left = stack_record["is_left"]

            if sample_indice is None:
                is_leaf = True
            else:
                n_samples = len(sample_indice)
                is_leaf = (depth >= self.max_depth or
                           n_samples < 2 * self.min_samples_leaf)

            if not is_leaf:
                split = self._node_split(sample_indice)
                if split is None:
                    is_leaf = True
                else:
                    impurity_improvement = impurity - split["impurity"]
                    is_leaf = (is_leaf or (impurity_improvement < (self.min_impurity_decrease + EPSILON)) or
                               (split["left"] is None) or (split["right"] is None))

            if is_leaf:
                node_id = self._add_node(parent_id, is_left, is_leaf, depth,
                                         None, None, impurity, sample_indice)
            else:
                node_id = self._add_node(parent_id, is_left, is_leaf, depth,
                                         split["feature"], split["threshold"], impurity, sample_indice)

                pending_node_list.append({"sample_indice": split["left"],
                                          "parent_id": node_id,
                                          "depth": depth + 1,
                                          "impurity": split["left_impurity"],
                                          "is_left": True})
                pending_node_list.append({"sample_indice": split["right"],
                                          "parent_id": node_id,
                                          "depth": depth + 1,
                                          "impurity": split["right_impurity"],
                                          "is_left": False})
        self.leaf_idx_list_ = [key for key, item in self.tree_.items() if item["is_leaf"]]
        return self

    def decision_rule(self, node_id):

        rule_dict = {}
        current_node = self.tree_[node_id]
        while True:
            if current_node["parent_id"] is None:
                break
            parent_node = self.tree_[current_node["parent_id"]]
            key = str(parent_node["feature"])
            if key not in rule_dict.keys():
                if current_node["is_left"]:
                    rule_dict.update({key: {"left": parent_node["threshold"]}})
                else:
                    rule_dict.update({key: {"right": parent_node["threshold"]}})
            else:
                if current_node["is_left"]:
                    if "left" not in rule_dict[key].keys():
                        rule_dict[key].update({"left": parent_node["threshold"]})
                    else:
                        rule_dict[key].update({"left": min(parent_node["threshold"], rule_dict[key]["left"])})
                else:
                    if "right" not in rule_dict[key].keys():
                        rule_dict[key].update({"right": parent_node["threshold"]})
                    else:
                        rule_dict[key].update({"right": max(parent_node["threshold"], rule_dict[key]["right"])})
            current_node = parent_node

        rule_list = []
        for key, item in rule_dict.items():
            rule = ""
            if "right" in item.keys():
                rule += str(round(item["right"], 3)) + "<"
            rule += "f" + key
            if "left" in item.keys():
                rule += "<=" + str(round(item["left"], 3))
            rule_list.append(rule)
        return rule_list

    def decision_path_indice(self, x, node_id):

        sample_indice = np.where(self.decision_path(x)[:, node_id - 1])[0]
        return sample_indice

    def decision_path(self, x):

        check_is_fitted(self)

        n_samples = x.shape[0]
        path_all = np.zeros((n_samples, self.node_count_))
        for node_id in self.leaf_idx_list_:
            path = []
            idx = node_id
            sample_indice = np.ones((x.shape[0],)).astype(bool)
            while True:
                path.append(idx - 1)
                current_node = self.tree_[idx]
                if current_node["parent_id"] is None:
                    break
                else:
                    parent_node = self.tree_[current_node["parent_id"]]
                    if current_node["is_left"]:
                        sample_indice = np.logical_and(sample_indice,
                                                       x[:, parent_node["feature"]] <= parent_node["threshold"])
                    else:
                        sample_indice = np.logical_and(sample_indice,
                                                       x[:, parent_node["feature"]] > parent_node["threshold"])
                idx = current_node["parent_id"]
            if sample_indice.sum() > 0:
                path_all[np.ix_(np.where(sample_indice)[0], path)] = 1
        return path_all

    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,),
        """
        n_samples = x.shape[0]
        pred = np.zeros((n_samples,))
        for node_id in self.leaf_idx_list_:
            idx = node_id
            sample_indice = np.ones((x.shape[0],)).astype(bool)
            while True:
                current_node = self.tree_[idx]
                if current_node["parent_id"] is None:
                    break
                else:
                    parent_node = self.tree_[current_node["parent_id"]]
                    if current_node["is_left"]:
                        sample_indice = np.logical_and(sample_indice,
                                                       x[:, parent_node["feature"]] <= parent_node["threshold"])
                    else:
                        sample_indice = np.logical_and(sample_indice,
                                                       x[:, parent_node["feature"]] > parent_node["threshold"])
                idx = current_node["parent_id"]
            if sample_indice.sum() > 0:
                predict_func = self.tree_[node_id]['predict_func']
                pred[sample_indice] = np.clip(np.dot(x[sample_indice, :], predict_func["coef"]),
                                              predict_func["min"], predict_func["max"]) + predict_func[
                                          "intercept"].ravel()
        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=False)
        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 MoGLMTreeRegressor(RegressorMixin, GLMTree, ModelBaseRegressor): """ A tree-based model that fits linear regression models in the leaves. This model recursively partitions the feature space and fits linear regression models in each leaf node. It combines the interpretability of decision trees with the flexibility of linear regression. Parameters ---------- name : str, default=None The name of the model. max_depth : int, default=3 The max number of depth. min_impurity_decrease : float, default=0 Minimum impurity decrease when splitting a node. min_samples_leaf : int, default=50 Minimum number of samples for constructing a leaf node. split_custom : dict, default=None The custom split points for each feature. n_screen_grid : int, default=1 The number of candidate split points for rough screening. n_feature_search : int, default=10 The number of candidate features selected by rough screening. n_split_grid : int, default=20 The number of candidate split points for fine-grained search. reg_lambda : float, default=0.1 The level of L1 regularization strength. clip_predict : bool, default=False Whether to clip the prediction results if it is outside the training data prediction. simplified : bool, default=True Whether to use partial linear regression for search the split feature and points. random_state : int, default=0 Determines random number generation for weights and bias initialization. Attributes ---------- tree_ : object The internal tree structure. """ def __init__(self, name: str = None, max_depth=3, min_samples_leaf=50, min_impurity_decrease=0, split_custom=None, n_screen_grid=1, n_feature_search=10, n_split_grid=20, clip_predict=False, reg_lambda=0.1, simplified=True, random_state=0): self.name = name super().__init__(base_estimator=xLinearRegression() if simplified else LinearRegression(), 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, 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"}} def _build_root(self): self.base_estimator.fit(self.x_, self.y_, self.sample_weight_) root_impurity = self._evaluate_estimator(self.base_estimator, self.x_, self.y_.ravel(), self.sample_weight_) return root_impurity def _build_leaf(self, sample_indice): EPSILON = 1e-7 mx = self.x_[sample_indice].mean(0) sx = self.x_[sample_indice].std(0) + EPSILON nx = (self.x_[sample_indice] - mx) / sx if self.reg_lambda_ is None: best_estimator = LinearRegression() elif isinstance(self.reg_lambda_, float): best_estimator = Lasso(alpha=self.reg_lambda_, precompute=False, random_state=self.random_state) elif isinstance(self.reg_lambda_, list) and len(self.reg_lambda_) > 1: best_estimator = LassoCV(alphas=self.reg_lambda_, cv=5, precompute=False, random_state=self.random_state) best_estimator.fit(nx, self.y_[sample_indice], self.sample_weight_[sample_indice]) best_estimator.coef_ = best_estimator.coef_ / sx best_estimator.intercept_ = best_estimator.intercept_ - np.dot(mx, best_estimator.coef_.T) if self.clip_predict: xmin = np.min(np.dot(self.x_[sample_indice], best_estimator.coef_)) xmax = np.max(np.dot(self.x_[sample_indice], best_estimator.coef_)) predict_func = {"coef": best_estimator.coef_, "intercept": best_estimator.intercept_, "min": xmin, "max": xmax} else: predict_func = {"coef": best_estimator.coef_, "intercept": best_estimator.intercept_, "min": -np.inf, "max": np.inf} best_impurity = self._get_loss(self.y_[sample_indice], best_estimator.predict(self.x_[sample_indice]), self.sample_weight_[sample_indice]) return predict_func, best_estimator, best_impurity @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 _evaluate_estimator(self, estimator, x, y, sample_weight=None): """method to calculate the MSE loss """ pred = estimator.predict(x) loss = self._get_loss(y, pred, sample_weight) return loss def _predict(self, X): """Predict function. Parameters ---------- X : np.ndarray Data features. Returns ------- np.ndarray The predicted label. """ X = check_array(X) check_is_fitted(self) return self.get_raw_output(X)
[docs] class MoGLMTreeClassifier(ClassifierMixin, GLMTree, ModelBaseClassifier): """ A tree-based model that fits logistic regression models in the leaves for binary classification. This model recursively partitions the feature space and fits logistic regression models in each leaf node. It combines the interpretability of decision trees with the flexibility of logistic regression. Parameters ---------- name : str, default=None Identifier name for the model instance. max_depth : int, default=3 Maximum depth of the tree. Controls model complexity. min_samples_leaf : int, default=50 Minimum number of 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 Dictionary mapping feature indices to custom split points. n_screen_grid : int, default=1 Number of grid points used in initial feature screening. n_feature_search : int, default=10 Number of top features to consider after screening. n_split_grid : int, default=20 Number of grid points to evaluate for splitting. reg_lambda : float, default=0.1 L1 regularization strength for leaf models. clip_predict : bool, default=False Whether to clip predictions to training data range. random_state : int, default=0 Random seed for reproducibility. Attributes ---------- tree_ : dict The fitted tree structure containing nodes and their parameters. leaf_estimators_ : dict Dictionary mapping leaf node IDs to their fitted logistic regression models. """ def __init__(self, name: str = None, max_depth=3, min_samples_leaf=50, min_impurity_decrease=0, split_custom=None, n_screen_grid=1, n_feature_search=10, n_split_grid=20, clip_predict=False, reg_lambda=0.1, random_state=0): self.name = name super().__init__(base_estimator=LogisticRegression(penalty=None, random_state=random_state), 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, random_state=random_state) def _more_tags(self): """ Internal function for skipping some sklearn estimator checks. """ return {"binary_only": True, "poor_score": True, "_xfail_checks": {"check_sample_weights_invariance": "zero sample_weight is not equivalent to removing samples"}} def _build_root(self): self.base_estimator.fit(self.x_, self.y_, self.sample_weight_) root_impurity = self._evaluate_estimator(self.base_estimator, self.x_, self.y_.ravel(), self.sample_weight_) return root_impurity def _build_leaf(self, sample_indice): EPSILON = 1e-7 if (self.y_[sample_indice].std() == 0) | (self.y_[sample_indice].sum() < 5) | ( (1 - self.y_[sample_indice]).sum() < 5): best_estimator = None predict_func = {"coef": np.zeros((self.x_.shape[1],)), "intercept": self.y_[sample_indice].mean(), "min": -np.inf, "max": np.inf} best_impurity = self._get_loss(self.y_[sample_indice], np.ones_like(self.y_[sample_indice]) * self.y_[sample_indice].mean(), self.sample_weight_[sample_indice]) else: # inverse the regularization as the "C" parameter in LogisticRegression is the inverse of regularization. if self.reg_lambda_ is None: best_estimator = LogisticRegression(C=1e7, random_state=self.random_state) elif isinstance(self.reg_lambda_, float): best_estimator = LogisticRegression(C=1 / self.reg_lambda_, random_state=self.random_state) elif isinstance(self.reg_lambda_, list) and len(self.reg_lambda_) > 1: best_estimator = LogisticRegressionCV(Cs=[1 / np.clip(reg, 1e-7, 1e7) for reg in self.reg_lambda_], penalty="l1", solver="liblinear", scoring="roc_auc", cv=5, random_state=self.random_state) mx = self.x_[sample_indice].mean(0) sx = self.x_[sample_indice].std(0) + EPSILON nx = (self.x_[sample_indice] - mx) / sx best_estimator.fit(nx, self.y_[sample_indice], self.sample_weight_[sample_indice]) best_estimator.coef_ = best_estimator.coef_ / sx best_estimator.intercept_ = best_estimator.intercept_ - np.dot(mx, best_estimator.coef_.T) xmin = np.min(np.dot(self.x_[sample_indice], best_estimator.coef_.ravel())) xmax = np.max(np.dot(self.x_[sample_indice], best_estimator.coef_.ravel())) if self.clip_predict: predict_func = {"coef": best_estimator.coef_.ravel(), "intercept": best_estimator.intercept_, "min": xmin, "max": xmax} else: predict_func = {"coef": best_estimator.coef_.ravel(), "intercept": best_estimator.intercept_, "min": -np.inf, "max": np.inf} best_impurity = self._get_loss(self.y_[sample_indice], best_estimator.predict_proba(self.x_[sample_indice])[:, 1], self.sample_weight_[sample_indice]) return predict_func, best_estimator, best_impurity @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 _evaluate_estimator(self, estimator, x, y, sample_weight=None): """method to calculate the cross entropy loss """ pred = estimator.predict_proba(x)[:, 1] loss = self._get_loss(y, pred, sample_weight) return loss 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