diff --git a/CHANGELOG.rst b/CHANGELOG.rst index ed5839252..31b940580 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -26,5 +26,6 @@ under development (0.3.2.dev0) - :bdg-danger:`Fix` Fixed unnecessary copy operations of X when only a slice view is needed (:gh:`646` by `Marc Hulcelle`_). - :bdg-secondary:`Maint` remove pinned poosh dependency (problem with somato dataset solved) (:gh:`670` by `Joseph Paillard`_). - :bdg-secondary:`Maint` remove extra term in variance of X-residual (DOCRT). See [Reid et al., A Study of Error Variance Estimation in Lasso Regression 2016](https://arxiv.org/pdf/1311.5274) for reference. (:gh:`649` by `Joseph Paillard`_). +- :bdg-success:`Feature` add leave-one-covariate-in (LOCI) method (:gh:`679` by `Marc Hulcelle`). - :bdg-danger:`Fix` fixed deprecated n_alphas with sklearn LassoCV, as well as deprecated penalty for LogisticRegressionCV (:gh:`690` by `Marc Hulcelle`) - :bdg-secondary:`Maint`: Fix conditional sampling test by verifying that sampler produces diverse samples. (:gh:`692` by `Joseph Paillard`_). diff --git a/docs/tools/references.bib b/docs/tools/references.bib index 8d342543f..8420824b4 100644 --- a/docs/tools/references.bib +++ b/docs/tools/references.bib @@ -483,3 +483,10 @@ @article{zhang2014confidence volume = {76}, year = {2014} } + +@inproceedings{qu2025tabicl, + title={Tab{ICL}: {A} Tabular Foundation Model for In-Context Learning on Large Data}, + author={Qu, Jingang and Holzm{\"u}ller, David and Varoquaux, Ga{\"e}l and Le Morvan, Marine}, + booktitle={International Conference on Machine Learning}, + year={2025} +} diff --git a/examples/plot_advanced_tabicl.py b/examples/plot_advanced_tabicl.py new file mode 100644 index 000000000..27b755ade --- /dev/null +++ b/examples/plot_advanced_tabicl.py @@ -0,0 +1,87 @@ +""" +Tabular Foundation Model TabICL +================================ + +In this example, we demonstrate how to use a tabular foundation model such as +TabICL [:footcite:t:`qu2025tabicl`] with a straightforward example on the California +Housing regression dataset. +""" + +# %% +# Loading data and TabICL +# ----------------------- +# We start by loading data from the California Housing dataset and fitting the +# TabICL model. If this is the first use, the model will be downloaded and cached +# for future use. The main advantage is that the model does not require +# hyperparameter tuning as a pre-trained transformer model. We recommend switching +# to GPU through the adequate class parameter for datasets that go over 10k +# samples depending on the number of features due to computation time increase. + +from sklearn.datasets import fetch_california_housing +from sklearn.model_selection import train_test_split +from sklearn.pipeline import make_pipeline +from sklearn.preprocessing import StandardScaler +from sklearn.utils import resample +from tabicl import TabICLRegressor + +dataset = fetch_california_housing() +X, y, feat_names = dataset.data, dataset.target, dataset.feature_names +X, y = resample( + X, y, replace=False, n_samples=1000, stratify=y, random_state=0 +) +X_train, X_test, y_train, y_test = train_test_split( + X, y, test_size=0.2, random_state=0 +) + +model = make_pipeline(StandardScaler(), TabICLRegressor(device="cpu")) +model.fit(X_train, y_train) +print(f"Model accuracy on the test data {model.score(X_test, y_test):.3}.") + +# %% +# For the moment, TabICL can only be used with PFI or CFI. TabICL only works with +# integer seeding, which is currently not supported by D0CRT. +# +# TabICL interfaces with HiDimStat the same way as any other Scikit-learn estimator, +# which lets us compute feature importance as easily as follows: + +import pandas as pd +from sklearn.metrics import mean_squared_error + +from hidimstat import CFI + +cfi = CFI(estimator=model, method="predict", loss=mean_squared_error) +cfi.fit(X_train, y_train) +importance = cfi.importance(X_test, y_test) +selection = cfi.fdr_selection(fdr=0.1) + +# %% +# Let's wrap the results in a Pandas dataframe and plot the importance +# and the selected variables. + +df = pd.DataFrame( + { + "feature": feat_names, + "importance": importance, + "selected": selection, + } +).sort_values("importance", ascending=False) + +import matplotlib.pyplot as plt +import seaborn as sns + +ax = sns.barplot( + data=df, + y="feature", + x="importance", + hue="selected", + palette="muted", + orient="h", +) +sns.despine() +plt.show() + +# %% +# As you can see from the timer under, even with 1000 samples and 9 features, +# running TabICL on the CPU is slow. We recommend to use it on larger datasets, +# be it in terms of samples and/or features and to run it on GPU for faster +# inferences. diff --git a/examples/plot_loci.py b/examples/plot_loci.py new file mode 100644 index 000000000..152277b34 --- /dev/null +++ b/examples/plot_loci.py @@ -0,0 +1,162 @@ +""" +Leave-One-Covariate-In (LOCI) feature importance with different regression models +================================================================================== + +This example demonstrates how to compare LOCI feature importance [:footcite:t:`Williamson_General_2023`] across different +predictive models on the same regression dataset. LOCI is model-agnostic and can be +applied to any predictive model. Here, we use a linear model, a random forest, a neural +network, and a support vector machine. We compare the models based on their predictive +performance (R2 score) and the LOCI feature importance they yield. +""" + +# %% +# Loading and preparing the data +# ------------------------------ +# We begin by simulating a regression dataset with 10 features, 5 of which +# are in the support set, meaning they contribute to generating the outcome. In this example, +# we use a simulated dataset to have access to the true support set of features and +# evaluate how well the different models identify these important features. +# The data is then split into training and test sets. These sets are used both to fit +# the predictive models and within the LOCI procedure, which refits models on subsets +# of features that exclude the feature of interest. + +from sklearn.datasets import make_regression +from sklearn.model_selection import train_test_split + +X, y, beta = make_regression( + n_samples=300, + n_features=10, + n_informative=5, + random_state=0, + coef=True, + noise=4.0, +) + +# We convert the coefficients of the data-generating process into a binary array +# indicating the true support set of features. +beta = beta != 0 + +X_train, X_test, y_train, y_test = train_test_split( + X, + y, + test_size=0.2, + random_state=0, + shuffle=True, +) + +# %% +# Fitting models and computing LOCI feature importance +# ---------------------------------------------------- +# We define a list of predictive models to compare. We use RidgeCV for linear +# regression, RandomForestRegressor for a tree-based model, MLPRegressor for a +# neural network, and SVR for a support vector machine, with RBF kernel. We then fit +# each model on the training data, compute the LOCI feature importance on the test +# data, and store the results in a DataFrame for comparison. +# Beware on the choice of hyperparameters of certain models to avoid over-fitting. + +import pandas as pd +from sklearn.ensemble import RandomForestRegressor +from sklearn.linear_model import RidgeCV +from sklearn.metrics import mean_squared_error +from sklearn.neural_network import MLPRegressor +from sklearn.svm import SVR + +from hidimstat import LOCI + +models_list = [ + RidgeCV(), + MLPRegressor( + hidden_layer_sizes=(8), + random_state=0, + max_iter=500, + learning_rate_init=0.1, + ), + SVR(kernel="linear"), + RandomForestRegressor( + n_estimators=50, max_depth=3, min_samples_leaf=10, random_state=0 + ), +] + +df_list = [] +for model in models_list: + # Fit the full model + model = model.fit(X_train, y_train) + loci = LOCI(model, method="predict", loss=mean_squared_error) + # Refit the model on a single feature / group of feature, and compute LOCI + # importance. This process is repeated for all features / groups of features to assess + # their individual contributions. + loci.fit(X_train, y_train) + importances = loci.importance(X_test, y_test) + df_list.append( + pd.DataFrame( + { + "feature": list(range(X.shape[1])), + "importance": importances, + "model": model.__class__.__name__, + "R2 score": model.score(X_test, y_test), + } + ) + ) + + +# %% +# The predictive performance of the models can be compared using their R2 scores. +# This helps assess how effectively each model captures the underlying data-generating +# process. + +df_plot = pd.concat(df_list) +df_plot.groupby("model").mean()["R2 score"].to_frame() + +# %% +# Visualization of LOCI feature importance +# ---------------------------------------- +# Finally, we visualize the LOCI feature importance for each model using a horizontal +# bar plot. The true support features are highlighted in the plot with a green shaded +# background. + +import matplotlib.pyplot as plt +import seaborn as sns + +ax = sns.barplot( + data=df_plot, + y="feature", + x="importance", + hue="model", + palette="muted", + orient="h", +) +sns.despine() + +for i, support in enumerate(beta): + if support != 0: + ax.axhspan( + i - 0.45, + i + 0.45, + color="tab:olive", + alpha=0.3, + zorder=-1, + label="True Support" if i == 1 else None, + ) +ax.legend() +plt.show() + +# %% +# The plot shows that the different models all identify the true support features and +# assign them higher importance scores. However, the magnitude of the importance scores +# varies across models. It can be observed that models with a greater predictive +# performance (higher R2 score) tend to assign higher importance scores to the true +# support features. + +# %% +# Conclusion +# ---------- +# LOCI is a simple and easily-interpretable method that can be used +# in parallel to LOCO to obtain a different interpretation of variable importance. +# Please keep in mind that models are fit on single features, which can lead to +# overfitting depending on the hyperparameters. One way to counteract this is to +# use LOCICV which can reduce noise for high-variance or unregularized models. + +# %% +# References +# ---------- +# .. footbibliography:: diff --git a/pyproject.toml b/pyproject.toml index 94bea6cb0..3a720f031 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -57,6 +57,7 @@ doc = [ "sphinx-gallery >= 0.17.0", "sphinx-prompt >= 1.9", "tqdm >= 4.66.3", + "tabicl >= 2.1.0", "mne >= 1.7.0", "sphinx-copybutton >= 0.5.2", "sphinx-design >= 0.5.0", diff --git a/src/hidimstat/__init__.py b/src/hidimstat/__init__.py index efc9ac217..4df065338 100644 --- a/src/hidimstat/__init__.py +++ b/src/hidimstat/__init__.py @@ -6,6 +6,7 @@ from .distilled_conditional_randomization_test import D0CRT, d0crt_importance from .ensemble_clustered_inference import CluDL, EnCluDL from .knockoffs import ModelXKnockoff, model_x_knockoff_importance +from .leave_one_covariate_in import LOCI, LOCICV, loci_importance from .leave_one_covariate_out import LOCO, LOCOCV, loco_importance from .permutation_feature_importance import PFI, PFICV, pfi_importance @@ -18,6 +19,8 @@ "CFI", "CFICV", "D0CRT", + "LOCI", + "LOCICV", "LOCO", "LOCOCV", "PFI", @@ -29,6 +32,7 @@ "cfi_importance", "d0crt_importance", "desparsified_lasso_importance", + "loci_importance", "loco_importance", "model_x_knockoff_importance", "pfi_importance", diff --git a/src/hidimstat/leave_one_covariate_in.py b/src/hidimstat/leave_one_covariate_in.py new file mode 100644 index 000000000..74781dcf8 --- /dev/null +++ b/src/hidimstat/leave_one_covariate_in.py @@ -0,0 +1,361 @@ +import numpy as np +import pandas as pd +from joblib import Parallel, delayed +from sklearn.base import check_is_fitted, clone, is_classifier, is_regressor +from sklearn.metrics import mean_squared_error + +from hidimstat._utils.docstring import _aggregate_docstring +from hidimstat._utils.utils import check_statistical_test +from hidimstat.base_perturbation import BasePerturbation, BasePerturbationCV + + +class LOCI(BasePerturbation): + """ + Leave-One-Covariate-In (LOCI) algorithm + + The model is re-fitted on each single feature/group of features. The importance is + then computed as the difference between the loss of an empty model (mean for regression, + and majority vote for classification) and the loss of the model on the single feature/group. + For more details, see :footcite:t:`ewald_2024`. + + Parameters + ---------- + estimator : sklearn compatible estimator + The estimator to use for the prediction. + method : str, default="predict" + The method to use for the prediction. This determines the predictions passed + to the loss function. Supported methods are "predict", "predict_proba" or + "decision_function". + loss : callable, default=mean_squared_error + The loss function to use when comparing the perturbed model to the full + model. + statistical_test : callable or str, default="ttest" + Statistical test function for computing p-values of importance scores. + features_groups: dict or None, default=None + A dictionary where the keys are the group names and the values are the + list of column names corresponding to each features group. If None, + the features_groups are identified based on the columns of X. + n_jobs : int, default=1 + The number of jobs to run in parallel. Parallelization is done over the + variables or groups of variables. + + References + ---------- + .. footbibliography:: + """ + + def __init__( + self, + estimator, + method: str = "predict", + loss: callable = mean_squared_error, + statistical_test="ttest", + features_groups=None, + n_jobs: int = 1, + ): + super().__init__( + estimator=estimator, + method=method, + loss=loss, + statistical_test=statistical_test, + features_groups=features_groups, + n_jobs=n_jobs, + ) + self._list_estimators = None + self._baseline_mean = None + + def fit(self, X, y): + """ + Fit a model for a single covariate/group of covariates. + + Parameters + ---------- + X : array-like of shape (n_samples, n_features) + The training input samples. + y : array-like of shape (n_samples,) + The target values. + + Returns + ------- + self : object + Returns the instance itself. + """ + super().fit(X, y) + # create a list of covariate estimators for each group if not provided + self._list_estimators = [ + clone(self.estimator) for _ in range(self.n_features_groups_) + ] + + # Parallelize the fitting of the covariate estimators + self._list_estimators = Parallel(n_jobs=self.n_jobs)( + delayed(self._joblib_fit_one_features_group)( + estimator, X, y, key_features_groups + ) + for key_features_groups, estimator in zip( + self.features_groups_.keys(), + self._list_estimators, + strict=False, + ) + ) + if self.method in ["predict_proba", "decision_function"]: + values, counts = np.unique(y, return_counts=True) + # Binary classification + if len(values) == 2: + values, counts = np.unique(y, return_counts=True) + self._baseline_mean = values[np.argmax(counts)] + # For multiclass classification, we take the marginal probability. + else: + self._baseline_mean = counts / y.shape[0] + elif self.method == "predict": + self._baseline_mean = np.mean(y) + return self + + def importance(self, X, y): + """ + Compute the importance scores for each group of covariates. + + Parameters + ---------- + X : array-like of shape (n_samples, n_features) + The input samples to compute importance scores for. + y : array-like of shape (n_samples,) + + Returns + ------- + importances_ : ndarray of shape (n_groups,) + The importance scores for each group of covariates. + A higher score indicates greater importance of that group. + + Attributes + ---------- + loss_reference_ : float + The loss of the model with the original (non-perturbed) data. + loss_ : dict + Dictionary with indices as keys and arrays of perturbed losses as values. + Contains the loss values for each permutation of each group. + importances_ : ndarray of shape (n_groups,) + The calculated importance scores for each group. + pvalues_ : ndarray of shape (n_groups,) + P-values from one-sided t-test testing if importance scores are + significantly greater than 0. + + Notes + ----- + The importance score for each group is calculated as the mean decrease in loss + when that feature group is included, compared to the null model. + A higher importance score indicates that including that feature group leads to + better model performance, suggesting those features are more important. + """ + self._check_fit() + self._check_compatibility(X) + statistical_test = check_statistical_test(self.statistical_test) + + if self.method in ["predict_proba", "decision_function"]: + values, _ = np.unique(y, return_counts=True) + # Binary classification + if len(values) == 2: + y_baseline = np.full_like(y, self._baseline_mean, dtype=int) + # Multiclass classification. + else: + y_baseline = np.full( + (y.shape[0], len(values)), self._baseline_mean + ) + elif self.method == "predict": + y_baseline = np.full_like(y, self._baseline_mean, dtype=float) + self.loss_reference_ = self.loss(y, y_baseline) + + y_pred = self._predict(X) + test_result = [] + self.loss_ = {} + for j, y_pred_j in enumerate(y_pred): + self.loss_[j] = np.array([self.loss(y, y_pred_j[0])]) + if np.all(np.equal(y.shape, y_pred_j[0].shape)): + test_result.append(y - y_pred_j[0]) + else: + test_result.append( + y - np.unique(y)[np.argmax(y_pred_j[0], axis=-1)] + ) + + self.importances_ = np.mean( + [ + self.loss_reference_ - self.loss_[j] + for j in range(self.n_features_groups_) + ], + axis=1, + ) + self.pvalues_ = statistical_test(np.array(test_result)).pvalue + assert self.pvalues_.shape[0] == y_pred.shape[0], ( + "The statistical test doesn't provide the correct dimension." + ) + return self.importances_ + + def _joblib_fit_one_features_group( + self, estimator, X, y, key_features_group + ): + """ + Fit the estimator on a group of covariates. + Used in parallel. + """ + if isinstance(X, pd.DataFrame): + X_j = X[self.features_groups_[key_features_group]] + else: + X_j = X[:, self.features_groups_[key_features_group]] + estimator.fit(X_j, y) + return estimator + + def _joblib_predict_one_features_group( + self, X, features_group_id, random_state=None + ): + """ + Predict the target feature for a single group of covariates. + Used in parallel. + """ + del random_state # not used (only there for API compatibility) + if isinstance(X, pd.DataFrame): + X_j = X[self.features_groups_[features_group_id]] + else: + # Since we don't have access to column names, we use the member _features_groups_ids + X_j = X[:, self._features_groups_ids[features_group_id]] + + y_pred_loci = getattr( + self._list_estimators[features_group_id], self.method + )(X_j) + + return [y_pred_loci] + + def _check_fit(self): + """Check that an estimator has been fitted after removing each group of + covariates. + """ + super()._check_fit() + check_is_fitted(self.estimator) + if self._list_estimators is None: + raise ValueError( + "The estimators require to be fit before to use them" + ) + for m in self._list_estimators: + check_is_fitted(m) + + +def loci_importance( + estimator, + X, + y, + method: str = "predict", + loss: callable = mean_squared_error, + features_groups=None, + test_statistic="ttest", + k_best=None, + percentile=None, + threshold_min=None, + threshold_max=None, + n_jobs: int = 1, +): + method = LOCI( + estimator=estimator, + method=method, + loss=loss, + statistical_test=test_statistic, + features_groups=features_groups, + n_jobs=n_jobs, + ) + method.fit_importance(X, y) + selection = method.importance_selection( + k_best=k_best, + percentile=percentile, + threshold_min=threshold_min, + threshold_max=threshold_max, + ) + return selection, method.importances_, method.pvalues_ + + +# use the docstring of the class for the function +loci_importance.__doc__ = _aggregate_docstring( + [ + LOCI.__doc__, + LOCI.__init__.__doc__, + LOCI.fit_importance.__doc__, + LOCI.importance_selection.__doc__, + ], + """ + Returns + ------- + selection : ndarray of shape (n_features,) + Boolean array indicating selected features (True = selected) + importances : ndarray of shape (n_features,) + Feature importance scores/test statistics. + pvalues : ndarray of shape (n_features,) + P-values computed for the marginal importance. + """, +) + + +class LOCICV(BasePerturbationCV): + """ + Leave-One-Covariate-IN (LOCI) algorithm with Cross-Validation. + + Parameters + ---------- + estimators: list of sklearn estimators or single sklearn estimator + Can be a list of fitted sklearn estimators (one per fold) or a single sklearn + estimator that will then be cloned and fitted on each fold. + cv: cross-validation generator + A cross-validation generator object (e.g., KFold, StratifiedKFold). + statistical_test : callable or str, default="nb-ttest" + Statistical test function to compute p-values from importance scores. + method : str, default="predict" + The method to use for the prediction. This determines the predictions passed + to the loss function. Supported methods are "predict", "predict_proba" or + "decision_function". + loss : callable, default=mean_squared_error + The loss function to use when comparing the perturbed model to the full + model. + features_groups: dict or None, default=None + A dictionary where the keys are the group names and the values are the + list of column names corresponding to each features group. If None, + the features_groups are identified based on the columns of X. + n_jobs : int, default=1 + The number of jobs to run in parallel. Parallelization is done over the folds. + + Attributes + ---------- + importance_estimators_ : list of LOCI instances + The LOCI instances fitted on each fold. + importances_ : ndarray of shape (n_groups, n_folds) + The calculated importance scores for each feature group and each fold. + Higher values indicate greater importance. + pvalues_ : ndarray of shape (n_groups,) + The p-values for the importance scores computed across folds. + estimators_ : list of sklearn estimators + List of fitted estimators for each fold. + test_train_frac_ : float + Fraction of test samples over train samples in each fold. Approximated as + 1 / (n_splits - 1). + """ + + def __init__( + self, + estimators, + cv, + statistical_test="nb-ttest", + method="predict", + loss=mean_squared_error, + features_groups=None, + n_jobs=1, + ): + super().__init__(estimators, cv, statistical_test, n_jobs) + self.method = method + self.loss = loss + self.features_groups = features_groups + + def _fit_single_split(self, estimator, X_train, y_train): + """Fit a LOCI instance on a single train/test split.""" + loci = LOCI( + estimator=estimator, + method=self.method, + loss=self.loss, + features_groups=self.features_groups, + n_jobs=1, # no parallelization inside the fold + ) + loci.fit(X_train, y_train) + return loci diff --git a/test/test_leave_one_covariate_in.py b/test/test_leave_one_covariate_in.py new file mode 100644 index 000000000..78bd39385 --- /dev/null +++ b/test/test_leave_one_covariate_in.py @@ -0,0 +1,253 @@ +from functools import partial + +import numpy as np +import pandas as pd +import pytest +from scipy.stats import ttest_1samp +from sklearn.datasets import make_classification +from sklearn.linear_model import LinearRegression, LogisticRegression, RidgeCV +from sklearn.metrics import log_loss +from sklearn.model_selection import KFold, train_test_split +from sklearn.preprocessing import OneHotEncoder + +from hidimstat import LOCI, LOCICV, loci_importance +from hidimstat._utils.scenario import multivariate_simulation +from hidimstat.base_perturbation import BasePerturbation +from hidimstat.statistical_tools.multiple_testing import fdp_power + + +def test_loci(): + """Test the Leave-One-Covariate-In algorithm on a linear scenario.""" + X, y, beta, _ = multivariate_simulation( + n_samples=100, + n_features=20, + support_size=4, + shuffle=False, + seed=42, + ) + important_features = np.where(beta != 0)[0] + non_important_features = np.where(beta == 0)[0] + + X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=0) + + regression_model = LinearRegression() + regression_model.fit(X_train, y_train) + + loci = LOCI( + estimator=regression_model, + method="predict", + features_groups=None, + n_jobs=1, + ) + + loci.fit( + X_train, + y_train, + ) + importance = loci.importance(X_test, y_test) + + assert importance.shape == (X.shape[1],) + assert ( + importance[important_features].mean() + > importance[non_important_features].mean() + ) + + # Same with groups and a pd.DataFrame + groups = { + "group_0": [f"col_{i}" for i in important_features], + "the_group_1": [f"col_{i}" for i in non_important_features], + } + X_df = pd.DataFrame(X, columns=[f"col_{i}" for i in range(X.shape[1])]) + X_train_df, X_test_df, y_train, y_test = train_test_split( + X_df, y, random_state=0 + ) + regression_model.fit(X_train_df, y_train) + loci = LOCI( + estimator=regression_model, + method="predict", + features_groups=groups, + n_jobs=1, + ) + loci.fit( + X_train_df, + y_train, + ) + importance = loci.importance(X_test_df, y_test) + + assert importance[0].mean() > importance[1].mean() + + # Classification case + y_clf = np.where(y > np.median(y), 1, 0) + _, _, y_train_clf, y_test_clf = train_test_split(X, y_clf, random_state=0) + logistic_model = LogisticRegression() + logistic_model.fit(X_train, y_train_clf) + + loci_clf = LOCI( + estimator=logistic_model, + method="predict_proba", + features_groups={ + "group_0": important_features, + "the_group_1": non_important_features, + }, + loss=log_loss, + n_jobs=1, + ) + loci_clf.fit( + X_train, + y_train_clf, + ) + importance_clf = loci_clf.importance(X_test, y_test_clf) + + assert importance_clf.shape == (2,) + assert importance_clf[0].mean() > importance_clf[1].mean() + + +def test_multiclass_loci(): + """Test LOCI in a multiclass classification setup.""" + important_features = 4 + n_features = 20 + + X, y = make_classification( + n_samples=100, + n_features=n_features, + n_informative=important_features, + n_classes=3, + random_state=42, + shuffle=False, + ) + groups = { + "important": list(range(important_features)), + "unimportant": list(range(important_features, n_features)), + } + + X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=0) + + logistic_model = LogisticRegression(random_state=0) + logistic_model.fit(X_train, y_train) + + loci_clf = LOCI( + estimator=logistic_model, + method="predict_proba", + features_groups=groups, + loss=log_loss, + n_jobs=1, + ) + loci_clf.fit( + X_train, + y_train, + ) + importance_clf = loci_clf.importance(X_test, y_test) + + assert importance_clf.shape == (2,) + assert importance_clf[0].mean() > importance_clf[1].mean() + + +def test_raises_value_error(): + """Test for error when model does not have predict_proba or predict.""" + X, y, _, _ = multivariate_simulation( + n_samples=150, + n_features=20, + support_size=4, + shuffle=False, + seed=42, + ) + + # Not fitted sub-model when calling importance and predict + with pytest.raises( + ValueError, match="This LOCI instance is not fitted yet" + ): + fitted_model = LinearRegression().fit(X, y) + loci = LOCI( + estimator=fitted_model, + method="predict", + ) + loci.importance(X, None) + + with pytest.raises( + ValueError, match="The estimators require to be fit before to use them" + ): + fitted_model = LinearRegression().fit(X, y) + loci = LOCI( + estimator=fitted_model, + method="predict", + ) + BasePerturbation.fit(loci, X, y) + loci.importance(X, y) + + with pytest.raises( + AssertionError, + match="The statistical test doesn't provide the correct dimension", + ): + fitted_model = LinearRegression().fit(X, y) + loci = LOCI( + estimator=fitted_model, + statistical_test=partial(ttest_1samp, popmean=0, axis=0), + ).fit(X, y) + loci.importance(X, y) + + +def test_loci_function(): + """Test the function of LOCI algorithm on a linear scenario.""" + X, y, beta, _ = multivariate_simulation( + n_samples=150, + n_features=20, + support_size=4, + shuffle=False, + seed=42, + ) + important_features = np.where(beta != 0)[0] + non_important_features = np.where(beta == 0)[0] + + X_train, _, y_train, _ = train_test_split(X, y, random_state=0) + + regression_model = LinearRegression() + regression_model.fit(X_train, y_train) + + _, importance, _ = loci_importance( + regression_model, + X, + y, + method="predict", + n_jobs=1, + ) + + assert importance.shape == (X.shape[1],) + assert ( + importance[important_features].mean() + > importance[non_important_features].mean() + ) + + +@pytest.mark.parametrize( + "n_samples, n_features, support_size, rho, seed, value, signal_noise_ratio, rho_serial", + [(500, 50, 5, 0.0, 0, 2.0, 8, 0.0)], + ids=["default data"], +) +def test_loci_cv(data_generator): + """ + Test that LOCI with cross-validated estimator works as expected. In particular, + - Empirical FDP is below the target FDR level + - Power is above 0.8, which is an arbitrary threshold + + Note: even though the only the expected FDP should be controlled, in practice + the simulation setting is simple enough to satisfy this stronger condition. + """ + X, y, important_features, _ = data_generator + + model = RidgeCV() + cv = KFold(n_splits=5, shuffle=True, random_state=0) + loci_cv = LOCICV( + estimators=model, + cv=cv, + n_jobs=5, + ) + loci_cv.fit(X, y) + loci_cv.importance(X, y) + + alpha = 0.2 + selected = loci_cv.fdr_selection(fdr=alpha) + gt_mask = np.zeros(X.shape[1], dtype=int) + gt_mask[important_features] = 1 + fdp, power = fdp_power(selected=selected, ground_truth=gt_mask) + assert fdp < alpha + assert power >= 0.8