---
title: "XAItest: Enhancing Feature Discovery with eXplainable AI"
shorttitle: XAItest
author: Ghislain FIEVET <ghislain.fievet@gmail.com>
package: XAItest
abstract: >
    XAItest is an R Package designed for integrating eXplainable Artificial
    Intelligence (XAI) techniques with traditional statistical analysis to
    enhance the discovery and interpretation of significant features within
    datasets. It applies p-values and feature importance metrics, such as
    SHAP, LIME or custom functions. A major XAItest tool is the generation
    of simulated data to establish significance thresholds, aiding in the
    direct comparison of feature importance with statistical p-values. The
    package aims at helping researchers and bioinformaticians seeking to
    uncover key features inside a dataset.
output:
    BiocStyle::html_document:
        toc: true
        toc_depth: 2
vignette: >
  %\VignetteIndexEntry{01_XAItest}
  %\VignetteEngine{knitr::knitr}
  %\VignetteEncoding{UTF-8}
bibliography: ../inst/REFERENCES.bib
---


## Introduction

In biological and biomedical research, the identification of key variables in
omics data—such as transcriptomes, proteomes, and methylomes—is generally
accomplished using statistical tests like the t-test or correlation analysis.
While these methods perform well for detecting classical mean differences,
they often fail to capture more complex scenarios, such as non-Gaussian data
distributions, interactions among multiple variables, or multimodal
distributions.

To address these challenges, the scientific literature proposes to use the
potential of explainable AI (XAI) to assist in the discovery of key features.
Examples of these techniques can be found in the analysis of methylation and
transcriptome data @Rajpal_2023, @Kumar_2023.

The XAItest package aims at proposing a versatile tool to run these types of
analysis easily. It calculates various p-value based statistics, such as the
t-test or correlation, depending on whether the target variable is categorical
or numerical. It also computes different XAI matrics such as feature
importance metrics, by default a built-in method from random forest, an
implementation of SHAP and of LIME, and the possibility to integrate custom
functions.

Establishing a significance level for a p-value is straightforward, often set
at 0.05 or 0.01. However, determining thresholds for the feature importance is
more complicated, and has to be computed at each experiment. Our package offers
an integrated solution to this challenge by generating simulated data designed
to meet selected p-value thresholds, thus setting significance thresholds for
the feature importance.

Lastly, our package offers an intuitive visualization of results through
color-coded tables, highlighting statistical metrics and their significance
for easy interpretation.


# Installation
```{r installation, eval=FALSE}
if (!require("BiocManager")) {
    install.packages("BiocManager")
}
BiocManager::install('XAItest')
```

# Generate benchmark datasets

The XAItest package provides several benchmark scenarios to evaluate the performance of statistical tests and machine learning methods. These scenarios include both classification and regression problems with different underlying patterns:

**Classification scenarios:**
- Scenario 1 : Variance Difference
- Scenario 2 : Bimodal Distribution
- Scenario 3 : XOR Interaction
- Scenario 4 : Concentric Circles

**Regression scenarios:**
- Scenario 5 : Parabolic Relationship
- Scenario 6 : Rising Sinusoid

These scenarios demonstrate cases where classical statistical tests (t-test, correlation) may fail to detect meaningful patterns, while machine learning approaches and XAI methods can successfully identify the underlying relationships.

```{r scenarios}
library(XAItest)
library(gridExtra)
library(ggplot2)

set.seed(12)

# Scenario 1: Variance Difference
df1 <- genScenario(1)
# Scenario 2: Bimodal Distribution
df2 <- genScenario(2)
# Scenario 3: XOR Interaction
df3 <- genScenario(3)
# Scenario 4: Concentric Circles
df4 <- genScenario(4)
# Scenario 5: Parabolic Relationship
df5 <- genScenario(5)
# Scenario 6: Rising Sinusoid
df6 <- genScenario(6)

# Create plots for each scenario
p1 <- ggplot(df1, aes(x = norm_noise01, y = diffVar01, color = y)) +
  geom_point() + scale_color_manual(values = c("orange", "purple")) +
  ggtitle("Scenario 1:\nVariance Difference")
p2 <- ggplot(df2, aes(x = norm_noise01, y = biDistrib, color = y)) +
  geom_point() + scale_color_manual(values = c("orange", "purple")) +
  ggtitle("Scenario 2:\nBimodal Distribution")
p3 <- ggplot(df3, aes(x = relVar01, y = relVar02, color = y)) +
  geom_point() + scale_color_manual(values = c("orange", "purple")) +
  ggtitle("Scenario 3:\nXOR Interaction")
p4 <- ggplot(df4, aes(x = relVar01, y = relVar02, color = y)) +
  geom_point() + scale_color_manual(values = c("orange", "purple")) +
  ggtitle("Scenario 4:\nConcentric Circles")
p5 <- ggplot(df5, aes(x = norm_noise01, y = y)) +
  geom_point() + ggtitle("Scenario 5:\nParabolic Relationship")
p6 <- ggplot(df6, aes(x = unif_noise01, y = y)) +
  geom_point() + ggtitle("Scenario 6:\nRising Sinusoid")

# Load required library for grid arrangement


# Arrange the plots in a 2x3 grid with width twice the height
grid.arrange(p1, p2, p3, p4, p5, p6, ncol = 3, nrow = 2,
             widths = rep(2, 3), heights = rep(1, 2))

```

# Add the `simThresh` simulated variable to the dataset

Objective of this section:

1. Add a simulated variable `simThresh` with XAItest::addSimThresh()
2. Verify that the p-value (t-test) is well calibrated
3. Train a decision tree and retrieve Gini importances
4. Use FI(simThresh) as threshold and identify features that exceed it

Brief explanation:
addSimThresh() adds a continuous variable "simThresh" constructed to
achieve a target p-value (default ~0.05) vs the target 'y'. It serves as a reference/benchmark for
comparing feature importances and p-values.

```{r addSimThresh}
set.seed(1)
# Add the simThresh variable to achieve a p-value = 0.001 between A and B
df_st <- addSimThresh(df1, pval_target = 0.001)
df_st[c(1:2,50:51),c(1:3,(ncol(df_st)-2):ncol(df_st))]
```

We verify that the p-value is indeed close to 0.001


```{r check pvalue}
grp_A <- df_st$simThresh[df_st$y == "A"]
grp_B <- df_st$simThresh[df_st$y == "B"]
tt <- t.test(grp_A, grp_B)
 tt$p.value
```

Now we train a decision tree model and compute Gini feature importances to identify which variables are most informative for classification. We'll use the `simThresh` variable's importance as a threshold to determine which features exceed this benchmark.

```{r build dt}
library(rpart)
fit <- rpart(y ~ ., data = df_st, method = "class")

imp_raw <- fit$variable.importance

# Significance threshold
fi_sim <- imp_raw[["simThresh"]]
fi_sim
```

Now that we have the threshold, we can identify which features exceed it.
```{r detect features}
# Feature Importance
imp_sorted <- sort(imp_raw, decreasing = TRUE)
imp_sorted
# Features above the threshold
vars_above <- setdiff(names(imp_sorted[imp_sorted > fi_sim]), "simThresh")
vars_above
```

These are indeed the key features from scenario 1: two classic mean difference variables (`diff_distrib01` and `diff_distrib02`), and the variance difference variable (`diffVar01`).

Now we will apply the same approach to a regression scenario using `df5`, which contains the parabolic relationship pattern.

```{r addSimThresh regress}
set.seed(1)
# Add the simThresh variable to achieve a p-value = 0.001 for the correlation with `y`
df_st <- addSimThresh(df5, pval_target = 0.001)
df_st[1:4,c(1:3,ncol(df_st))]

tt <- cor.test(df_st$y, df_st$simThresh )
tt$p.value

fit <- rpart(y ~ ., data = df_st)

imp_raw <- fit$variable.importance
imp_raw

# Significance threshold
fi_sim <- imp_raw[["simThresh"]]
fi_sim

# Feature Importance
imp_sorted <- sort(imp_raw, decreasing = TRUE)
imp_sorted
# Features above the threshold
vars_above <- setdiff(names(imp_sorted[imp_sorted > fi_sim]), "simThresh")
vars_above
```

The feature `norm_noise01` is indeed the feature used to build the parabolic `y` feature.


# Get feature importance of a SummarizedExperiment dataset

```{r load libs, message=F}
# Load the libraries
library(ggforce)
library(SummarizedExperiment)
```
We consider the RNA-seq data from the *airway* dataset, which contains read
counts per gene for airway smooth muscles stored in a
*RangedSummarizedExperiment* object. The metadata contains a column *dex*,
which indicates whether the subjects are treated or untreated. The *XAI.test*
function, by default, builds random forest models to predict a target, here
the 'dex' treatment status, and calculates the associated gene feature
importances.

```{r}
data(airway, package="airway")
se <- airway
```

```{r first XAI.test, warning=FALSE}
results <- XAI.test(se[1:100,], y = "dex", verbose = TRUE, simData = TRUE)
```

```{r}
results
```

# Detect new key features on a classification dataset
On a toy dataset we will investigate which features can be identified using a
t-test or linear modeling, and which can be detected through feature importance
metrics.

## Load simulated classification dataset
The dataset contains 50 samples of **A**, and 50 samples of **B**, with 10
noise features et 3 key features:
- **norm_noise01**, **norm_noise02** ... 10 features of normal distribution
    noise
- **diff_distrib01** different normal distributions between "A" and "B"
- **diff_distrib02** different normal distributions between "A" and "B"
- **bidistrib01** normal bidistribution for "A" and normal distribution for "B"

The t-test can detect **diff_distrib01** and **diff_distrib02** but not
**bidistrib01**. 

```{r load classif simu data}
se_path <- system.file("extdata", "seClassif.rds", package="XAItest")
dataset_classif <- readRDS(se_path)

data_matrix <- assay(dataset_classif, "counts")
data_matrix <- t(data_matrix)
metadata <- as.data.frame(colData(dataset_classif))
df_simu_classif <- as.data.frame(cbind(data_matrix, y = metadata[['y']]))
for (col in names(df_simu_classif)) {
    if (col != 'y') {
        df_simu_classif[[col]] <- as.numeric(df_simu_classif[[col]])
    }
}
```

Example with a dataframe.

```{r load dataframe ex classif simu data}
df_path <- system.file("extdata", "dfClassif.txt", package="XAItest")
dataset_classif_df <- read.table(df_path)
df_simu_classif_df <- dataset_classif
```

We can see below how class **A** and class **B** are distributed accross the
features: **norm_noise01**, **norm_noise02**, **norm_noise03**,
**norm_noise04**, **diff_distrib01**, **diff_distrib02** and **biDistrib**.

```{r fig.width=12, fig.height=5}
p1 <- ggplot(df_simu_classif, aes(x=norm_noise01, y=norm_noise02,
            color=y)) +
    geom_point() +
    ggtitle("Noise features\nnorm_noise01 vs norm_noise02") +
    theme_bw()
p2 <- ggplot(df_simu_classif, aes(x=diff_distrib01, y=diff_distrib02,
            color=y)) +
    geom_point() +
    ggtitle("Normal distributions\ndiff_distrib01 vs diff_distrib02") +
    theme_bw()
p3 <- ggplot(df_simu_classif, aes(x=norm_noise03, y=biDistrib, color=y)) +
    geom_point() +
    ggtitle("Normal bidistribution\nbidistrib01 vs norm_noise03") +
    theme_bw()
grid.arrange(p1, p2, p3, ncol = 3, nrow = 1)
```

We employ the XAI.test function to evaluate the p-values and feature importance
of each feature. The parameter *simData* is set to TRUE, which instructs the
function to generate a new feature named *simThresh*, designed to achieve a
specified p-value (column *ttest_adjPval*), set at *simPvalTarget = 0.01*.
The purpose of *simThresh* is to establish a significance threshold for assessing
feature importance values.

```{r XAI.test 1 , warning=FALSE}
set.seed(123)
objXAI <- XAI.test(dataset_classif, "y", simData = TRUE, simPvalTarget = 0.01)
```

```{r}
head(getMetricsTable(objXAI))
```

The *mapPvalImportance* function reveals the significance of the feature
importance by comparison with the p-values.

Display as a data.frame:
```{r mapPvalImportance1}
mpi <- mapPvalImportance(objXAI)
mpi$df
```

Display as a datatable:
```{r mapPvalImportance5, eval=FALSE}
mpi$dt
```

![ ](im/dt1.png)

We can see that the **biDistrib** feature **is not detected by the p-values**,
but is detected by the **feature importance** metrics.

## Plot each model

The *XAI.test* function builds several predictive models.
In order to provide a p-value, linear regression (lm) builds a model that we
can access. Similarly, to give feature importances, methods like the random
forest, SHAP, and LIME also build models that are accessible

The *plotModel* function is used to visualize the results of these predictive
models. This visualization helps determine which classes are well-predicted by
which model. This function aids in more effectively interpreting the model's
performance in each situation.

In this section, we aim to evaluate the performance of each model after
removing the *diff_distrib01* and *diff_distrib02* features from the dataset,
leaving only noise and the *biDistrib* features.

```{r xai2, warning=FALSE}
set.seed(123)
objXAI <- XAI.test(df_simu_classif[,setdiff(colnames(df_simu_classif),
                                    c("diff_distrib01", "diff_distrib02"))],
                  simData=TRUE)
head(getMetricsTable(objXAI))
```

We can see below that the linear model predictions contain many errors, while
the random forest, along with the models used for SHAP and LIME, are much more
precise.

```{r fig.width=10, fig.height=10}
p1 <- plotModel(objXAI, "lm_pval", "biDistrib", "simThresh")
p2 <- plotModel(objXAI, "RF_feat_imp", "biDistrib", "simThresh")
p3 <- plotModel(objXAI, "SHAP_feat_imp", "biDistrib", "simThresh")
p4 <- plotModel(objXAI, "LIME_feat_imp", "biDistrib", "simThresh")
grid.arrange(p1, p2, p3, p4, ncol = 2, nrow = 2)
```

Display as a data.frame:
```{r mapPvalImportance2}
mpi <- mapPvalImportance(objXAI)
head(mpi$df)
```

Display as a datatable:
```{r mapPvalImportance15, eval=FALSE}
mpi$dt
```
![ ](im/dt2.png)

# Detect new key features on a regression dataset

First, we load a regression toy example of 100 samples, with 5 noise features,
and a *y* feature computed from one of the noise feature. We will investigate
if the feature used to compute *y* can be identified using correlation, linear
modeling, or feature importance metrics.

## Load regression simulated dataset
- **norm_noise1**, **norm_noise2**: normal distribution noise features
- **unif_noise1**, **unif_noise2** and **unif_noise3**: uniform distribution
    noise features
- **y** feature is aparabole computed from the **norm_noise2** values.

### Function of parabole
This function is invisible to the correlation test and is used to create the
output `df_simu_regr$y <- transfo_parab(df_simu_regr$norm_noise2)`.

```{r}
transfo_parab <- function(xs){
    x1 <- min(xs)
    x2 <- max(xs)
    h <- (x1 + x2) / 2
    k <- max(xs)/2
    a <- k / ((h - x1) * (h - x2))
    
    y <- a * (xs - x1) * (xs - x2)
    
    return (y)
}

```


```{r load regress simu data}
se_path <- system.file("extdata", "seRegress.rds", package="XAItest")
dataset_regress <- readRDS(se_path)

data_matrix <- assay(dataset_regress, "counts")
data_matrix <- t(data_matrix)
metadata <- as.data.frame(colData(dataset_regress))
df_simu_regr <- as.data.frame(cbind(data_matrix, y = metadata[['y']]))
for (col in names(df_simu_regr)) {
    if (col != 'y') {
        df_simu_regr[[col]] <- as.numeric(df_simu_regr[[col]])
    }
}
```

Example with a dataframe.
```{r load dataframe regress simu data}
df_path <- system.file("extdata", "dfRegress.txt", package="XAItest")
dataset_regress_df <- read.table(df_path)
df_simu_regr_df <- dataset_regress
```

```{r}
ggplot(df_simu_regr, aes(x=norm_noise2, y=y)) + geom_point() + theme_bw()
```

We employ the *XAI.test* function to evaluate the p-values and feature
importance of each feature. The parameter *simData* is set to TRUE, which
instructs the function to generate a new feature named *simThresh*, which is
designed to achieve a specified p-value (column *ttest_adjPval*), set at
*simPvalTarget = 0.01*. The purpose of *simThresh* is to establish a significance
threshold for assessing feature importance values.

```{r xai3, warning=FALSE}
set.seed(123)
regr_results <- XAI.test(dataset_regress, "y",
                            simData=TRUE, simPvalTarget = 0.01)
getMetricsTable(regr_results)
```

Display as a data.frame:
```{r mapPvalImportance3}
mpi <- mapPvalImportance(regr_results, refPvalColumn = "cor_adjPval", refPval = 0.01)
head(mpi$df)
```
Display as a datatable:
```{r mapPvalImportance25, eval=FALSE}
mpi$dt
```
![ ](im/dt3.png)

The relationship between the *norm_noise2* and *y* features is not apparent
through correlation and linear modeling statistical tests, but is revealed by
the **Random Forest built-in feature importance**, as well as its **SHAP** and
**LIME** values.

## Plot each model
Similar to the classification case, the use of the *XAI.test* function creates
predictive models that we can visualize with the *plotModel* function. This
provides insights into why the *norm_noise2* feature was detected by feature
importance metrics and not by linear modeling.

```{r}
regr_results@args$modelType
```

```{r fig.width=10, fig.height=5}
p1 <- plotModel(regr_results, "lm_pval", "norm_noise2")
p2 <- plotModel(regr_results, "SHAP_feat_imp", "norm_noise2")
grid.arrange(p1, p2, ncol = 2, nrow = 1)
```

Use the *modelsOverview* function to quickly review the performance of the
models.

```{r}
modelsOverview(regr_results)
```

# Using custom functions

Here is `fctParamFIPV()`, an example of a custom function used in the preprint computations.
This function still requires further improvements before being production-ready.
Set `simData=T` to add the `simThresh` variable.

```{r cust fct}
script_path <- system.file("fctParamFIPV.R", package = "XAItest")
source(script_path)

xai_res <- XAI.test(df1, y="y", simData=F, customFIPV=list("pimp"=function(data, ...) {
    fctParamFIPV(data=data, fi_algo="gini", pv_algo="pimp", ml_algo="dt",
    fi_transfo="none", n_perm=20, shap_nsim = 10, verbose=0, ...)
}), defaultMethods=c())

xai_res
```

The number of permutations in PIMP and mProbes is controlled by `n_perm`, at least 100 recommended.
The number of probes in mProbes, or the number of combinations for SHAP and LIME, is specified by `shap_nsim`.

You can see the compute times for each column:
```{r times}
xai_res@computeTimes
```


The available arguments combinations are:

- `ml_algo="dt", fi_algo="gini", pv_algo="pimp"`
- `ml_algo="dt", fi_algo="gini", pv_algo="mprobes"`
- `ml_algo="dt", fi_algo="shap", pv_algo="pimp"`
- `ml_algo="dt", fi_algo="shap", pv_algo="mprobes"`
- `ml_algo="dt", fi_algo="lime", pv_algo="pimp"`
- `ml_algo="dt", fi_algo="lime", pv_algo="mprobes"`
- `ml_algo="rf", fi_algo="gini", pv_algo="pimp"`
- `ml_algo="rf", fi_algo="gini", pv_algo="mprobes"`
- `ml_algo="rf", fi_algo="accuracy", pv_algo="pimp"`
- `ml_algo="rf", fi_algo="accuracy", pv_algo="mprobes"`
- `ml_algo="rf", fi_algo="shap", pv_algo="pimp"`
- `ml_algo="rf", fi_algo="shap", pv_algo="mprobes"`
- `ml_algo="rf", fi_algo="lime", pv_algo="pimp"`
- `ml_algo="rf", fi_algo="lime", pv_algo="mprobes"`
- `ml_algo="svm+vallinadot", fi_algo="shap", pv_algo="pimp"`
- `ml_algo="svm+vallinadot", fi_algo="shap", pv_algo="mprobes"`
- `ml_algo="svm+vallinadot", fi_algo="lime", pv_algo="pimp"`
- `ml_algo="svm+vallinadot", fi_algo="lime", pv_algo="mprobes"`
- `ml_algo="svm-rbfdot", fi_algo="shap", pv_algo="pimp"`
- `ml_algo="svm-rbfdot", fi_algo="shap", pv_algo="mprobes"`
- `ml_algo="svm-rbfdot", fi_algo="lime", pv_algo="pimp"`
- `ml_algo="svm-rbfdot", fi_algo="lime", pv_algo="mprobes"`
- `ml_algo="mlp", fi_algo="olden", pv_algo="pimp"`
- `ml_algo="mlp", fi_algo="olden", pv_algo="mprobes"`
- `ml_algo="mlp", fi_algo="shap", pv_algo="pimp"`
- `ml_algo="mlp", fi_algo="shap", pv_algo="mprobes"`
- `ml_algo="mlp", fi_algo="lime", pv_algo="pimp"`
- `ml_algo="mlp", fi_algo="lime", pv_algo="mprobes"`



```{r}
sessionInfo()
```

# References