---
title: "findIPs: Detect Influential Points for Feature Rankings"
author:
- name: Shuo Wang
  affiliation: 
  - Medical Faculty Heidelberg, Heidelberg University
  - Institute of Medical Biometry and Statistics, University of Freiburg
  email: wangsures@foxmail.com
- name: Junyan Lu
  affiliation: Medical Faculty Heidelberg, Heidelberg University
date: "`r format(Sys.time(), '%B %d, %Y')`"

output:
   BiocStyle::html_document: default
   BiocStyle::pdf_document: default
link-citations: yes
header-includes:
    - \usepackage{setspace}
    - \doublespacing
    
vignette: >
  %\VignetteIndexEntry{Introduction to package findIPs}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

# Introduction

Feature rankings are important in analyzing high-throughput data, particularly 
for bioinformatics studies. These rankings are typically obtained by calculating
the marginal correlations with the outcomes. The ordered feature list reflects 
the importance of features, which plays a vital role in guiding subsequent 
research. For instance, researchers usually focus on a small subset of important
features that are associated with research objectives. However, the feature 
ranking can be distorted by a single case. The case exerts abnormal influence 
on the feature ranking is termed as influential points (IPs). The presence of 
IPs renders the feature ranking unreliable, consequently affecting the 
subsequent analysis based on feature rankings.

The *findIPs* R package is specifically designed to detect IPs for feature 
rankings. The method utilized in this package is based on the leave-one-out 
strategy, which involves comparing the rank difference between the original 
ranking and a new ranking that is obtained by removing a single observation (1). 
The  new rankings are leave-one-out rankings. The rank difference obtained 
through this comparison helps to determine the influence of the deleted 
observation.

The whole process can be divided into three steps,  

Step 1, generate the original and leave-one-out rankings using a feature 
ranking method, such as t-test. A dataset with n cases will result in one 
original ranking and n leave-one-out feature rankings.

Step 2, calculate rank changes. It is advisable to use top-prioritized 
weights when comparing ranks.

Step 3, calculate the cumulative rank changes for each observation. A diagnostic
check is also required to identify any potential influential points. 


# Installation

The *findIPs* package can be installed from  Bioconductor using the 
following commands:

```{r, eval = FALSE}

if (!requireNamespace("BiocManager", quietly = TRUE))
     install.packages("BiocManager")
BiocManager::install("findIPs")

```

# Dataset

The *findIPs* package includes the miller05 microarray data (2). The data 
contains 236 samples and 1,000 genes. It has two types of responses: binary 
and survival. The binary response is classified based on the p53 mutation: 58 
cases with a p53 mutant and 193 cases with the wild-type p53 mutation. The 
survival response is associated with a total of 55 recorded events. 

```{r}
library(findIPs)
data(miller05)
X <- miller05$X
y <- miller05$y
surv <- miller05$surv
```

# Detect IPs using *getdrop1ranks()* and *sumRanks()*  

We use a simple example where features are ranked based on t.test to demonstrate 
the use of *findIPs* package for IPs detection. We use *getdrop1ranks()* to 
derive the original ranking and leave-one-out rankings. Features are simply 
ranked according to the p.values of t.test. Of note, the rank criteria is the 
p.value if *fun = "t.test"*. P.values are ranked in ascending order by 
specifying *decreasing = FALSE*. We select the top 100 important features in 
the original ranking. The function returns an object containing a vector of 
original ranking (origRank) and a matrix of leave-one-out rankings (drop1Rank).

```{r}

obj <- getdrop1ranks(X, y,
                     fun = "t.test",
                     decreasing = FALSE,
                     topN = 100)
str(obj)

```

After obtaining the original ranking and leave-one-out rankings using the 
*getdrop1ranks()* function, we use the *sumRanks()* function to compute the 
distance between them. This function provides three methods for comparing ranks: 
unweighted, weighted Spearman, and method with adaptive weights. The unweighted 
method directly compares the ranks and assumes that all ranks have equal 
importance. However, this is not always the case as the top-ranked methods are 
usually more important. The weighted Spearman and adaptive weights methods 
address this issue by emphasizing the importance of the top-ranked methods (3). 
The adaptive weights method can further adjust the weights based on the 
distribution of rank changes in the data.


```{r}

results <- sumRanks(origRank = obj$origRank,
                    drop1Rank = obj$drop1Rank,
                    topN = 100,
                    method = "adaptive")

str(results)

```

The outputs of *sumRanks()* are diverse across the selected methods. For 
*method = "adaptive"*, *sumRanks()* returns a list with the following elements:

1, kappa, the shape parameter of the adaptive weights method;  
2, score, the accumulated weighted rank changes, reflecting the influence of 
each sample;  
3, origRank, the original ranking;  
4, drop1Rank, the leave-one-out rankings;   
5, origRankWeighted, weighted original ranking;  
6, drop1RankWeighted, weighted leave-one-out rankings.  

However, if the method is "weightedSpearman" or "unweighted", the function will 
only return three elements: "score", "origRank", and "drop1RankWeighted". The 
elements "kappa", "origRankWeighted", and "drop1RankWeighted" will not be 
returned.

# Use *findIPs()* to detect IPs in one-step

*findIPs()* combines *getdrop1ranks()* and *sumRanks()* into one step. The 
output is identical to that using the two-step process.

```{r}
results.ipsr <- findIPs(X, y, 
                        fun = "t.test",
                        decreasing = FALSE,
                        method = "adaptive")

identical(results, results.ipsr)
```
# Results visualization 

*findIPs* package offers three visualization functions: *plotIPs()*, 
*plotRankScatters()*, and *plotAdaptiveWeights()*. *plotIPs()* can directly 
utilize the output of *findIPs()* or *sumRanks()* to create a lollipop plot that 
displays the influence of each case. In Figure 1, we can see that the 
observation 68 (obs68) seems to be more influential on the final results. 
However, the difference between obs68 and the other observations is not that 
distinct, indicating a lower possibility of the presence of an influential 
observation.

```{r, fig.width = 7, fig.height = 4, fig.cap = "Figure 1, the influence scores for each observation"}
par(mar = c(4, 4, 2, 2))
plotIPs(results.ipsr, topn = 5, ylim = c(0, 8))

```

In addition to the lollipop, *findIPs* also provides a simple visualization 
function *plotRankScatters()* that exhibits the distribution of rank changes 
using a scatter plot (Figure 2). Alike *plotIPs()*, *plotRankScatters()* simply 
uses the output of *findIPs()* or *sumRanks()*. According to Figure 2, we can 
observe more rank changes in the tail side, but less changes in the head. The 
black points denote the rank changes caused by the most influential case.  

```{r, fig.width = 7, fig.height = 4, fig.cap = "Figure 2, the distribution of rank changes"}
par(mar = c(4, 4, 2, 2))
plotRankScatters(results.ipsr)

```

The *plotAdaptiveWeights()* function aims to visualize weight function if 
adaptive weights are used for rank comparison, that is *method = "adaptive"* 
for *findIPs()* or *sumRanks()*. The argument *kappa* refers to the shape 
parameter of the weight function. Here, the optimized kappa is 0.023 (Figure 3). 
n is the length of the feature list. We select the top 100 features, hence, 
*n = 100*. We can observe that more weights are allocated to the top-ranked 
features.  

```{r, fig.width = 7, fig.height = 4, fig.cap = "Figure 3, the weight function of the adaptive weights"}
par(mar = c(4, 4, 2, 2))
plotAdaptiveWeights(results.ipsr$kappa, n = 100, type = "line")

```

# Use findIPs in survival data

For survival analysis, we offer the option to rank features using univariate 
Cox model by setting *fun = "cox"*. The features are ranked in ascending order 
based on their P-values.

```{r, fig.width = 7, fig.height = 4, fig.cap = "Figure 4, IPs detection for survival data"}
par(mar = c(4, 4, 2, 2))
results.cox <- findIPs(X, surv, 
                       fun = "cox",
                       decreasing = FALSE,
                       method = "adaptive")

plotIPs(results.cox)
```


# Customize the rank criteria

In addition to the provided ranking criteria in *findIPs()* or *sumRanks()*, 
which includes "t.test", "cox", "log2fc", and "kruskal.test". We can also rank 
features based on a specified rank criterion. To this end, we can pass a function 
to the *fun* argument.  The function should take *x* and *y* as inputs and 
output the rank criterion, such as p-values.

As an example, we can rank features based on the p-values obtained from the 
kruskal.test. We can either specify *fun = "kruskal.test"*, as this test has 
been implemented in the package, or define our own function passing to 
*getdrop1ranks*. Both methods produce the same results.

```{r}
fun <- function(x, y){
  kruskal.test(x, y)$p.value
}

kruskal.test1 <- getdrop1ranks(X, y, 
                               fun = fun, 
                               decreasing = FALSE)

kruskal.test2 <- getdrop1ranks(X, y, 
                               fun = "kruskal.test", 
                               decreasing = FALSE)

identical(kruskal.test1, kruskal.test2)

```

# The choice of rank comparison methods

*findIPs* provides three rank comparison methods: unweighted, weighted Spearman, 
and adaptive weights. We recommend using the adaptive weights. Here, we compare 
the three methods. Weighted Spearman and the adaptive weights method 
demonstrates similar results. But both are different from the unweighted method. 

```{r, fig.width = 9, fig.height = 3, fig.cap = "Figure 5, IPs detection using three rank comparison methods"}
re.unweighted <- findIPs(X, y, 
                         fun = "t.test",
                         decreasing = FALSE,
                         method = "unweighted")
re.weighted <- findIPs(X, y, 
                       fun = "t.test",
                       decreasing = FALSE,
                       method = "weightedSpearman")
re.adaptive <- findIPs(X, y, 
                       fun = "t.test",
                       decreasing = FALSE,
                       method = "adaptive")

par(mfrow = c(1, 3), mar = c(4, 4, 2, 2))
plotIPs(re.unweighted)
mtext("A, unweighted", side = 3, line = 0, adj = 0)
plotIPs(re.weighted)
mtext("B, weighted Spearman", side = 3, line = 0, adj = 0)
plotIPs(re.adaptive)
mtext("C, adaptive weights", side = 3, line = 0, adj = 0)

```

# References {.unnumbered}

1, Wang, Shuo, and Junyan Lu. "Detect influential points of feature rankings." 
arXiv preprint arXiv:2303.10516 (2023).  
2, Miller, Lance D., et al. "An expression signature for p53 status in human 
breast cancer predicts mutation status, transcriptional effects, and patient 
survival." Proceedings of the National Academy of Sciences 102.38 (2005): 
13550-13555.doi:10.1073pnas.0506230102  
3, Da Pinto Costa, Joaquim; Soares, Carlos (2005): A weighted rank measure of 
correlation. In Australian & New Zealand Journal of Statistics 47 (4), pp. 
515–529.  


# Session info {.unnumbered}

```{r, echo = F}
sessionInfo()
```


