---
title: "ROCaggregator use case"
author: "Pedro Mateus"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{ROCaggregator use case}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

## Introduction

The ROCaggregator package allows you to aggregate multiple ROC (Receiver 
Operating Characteristic) curves. One of the scenarios where it can be helpful 
is in federated learning. Evaluating a model using the ROC AUC (Area Under the 
Curve) in a federated  learning scenario will require to evaluating the model 
against data from different sites. This will eventually lead to partial ROCs 
from each site which can be aggregated to obtain a global metric to evaluate 
the model.

```{r intro}
library(ROCaggregator)
```

## Set-up

For this use case, we'll be using some external packages:
- the `ROCR` to compute the ROC at each node and to validate the AUC obtained;
- the `pracma` package to compute the AUC using the trapezoidal method;
- the `stats` package to create a linear model;

The use case will consist of 3 nodes with horizontally partitioned data. A 
linear model will be trained with part of the data and tested at each node, 
generating a ROC curve for each. 

To compute the aggregated ROC, each node will have to provide:
- the ROC (consisting of the false positive rate and true positive rate);
- the thresholds/cutoffs used (in the same order as the ROC);
- the total number of negative labels in the dataset;
- the total number of samples in the dataset;

```{r setup}
library(ROCR)
library(pracma)
library(stats)

set.seed(13)

create_dataset <- function(n){
  positive_labels <- n %/% 2
  negative_labels <- n - positive_labels

  y = c(rep(0, negative_labels), rep(1, positive_labels))
  x1 = rnorm(n, 10, sd = 1)
  x2 = c(rnorm(positive_labels, 2.5, sd = 2), rnorm(negative_labels, 2, sd = 2))
  x3 = y * 0.3 + rnorm(n, 0.2, sd = 0.3)
  
  data.frame(x1, x2, x3, y)[sample(n, n), ]
}

# Create the dataset for each node
node_1 <- create_dataset(sample(300:400, 1))
node_2 <- create_dataset(sample(300:400, 1))
node_3 <- create_dataset(sample(300:400, 1))

# Train a linear model on a subset
glm.fit <- glm(
  y ~ x1 + x2 + x3,
  data = rbind(node_1, node_2),
  family = binomial,
)

get_roc <- function(dataset){
  glm.probs <- predict(glm.fit,
                       newdata = dataset,
                       type = "response")
  pred <- prediction(glm.probs, c(dataset$y))
  perf <- performance(pred, "tpr", "fpr")
  perf_p_r <- performance(pred, "prec", "rec")
  list(
    "fpr" = perf@x.values[[1]],
    "tpr" = perf@y.values[[1]],
    "prec" = perf_p_r@y.values[[1]],
    "thresholds" = perf@alpha.values[[1]],
    "negative_count"= sum(dataset$y == 0),
    "total_count" = nrow(dataset),
    "auc" = performance(pred, measure = "auc")
  )
}

# Predict and compute the ROC for each node
roc_node_1 <- get_roc(node_1)
roc_node_2 <- get_roc(node_2)
roc_node_3 <- get_roc(node_3)
```

## Aggregating the ROC from each node

Obtaining the required inputs from each node will allow us to compute the 
aggregated ROC and the corresponding AUC.

```{r aggregating}
# Preparing the input
fpr <- list(roc_node_1$fpr, roc_node_2$fpr, roc_node_3$fpr)
tpr <- list(roc_node_1$tpr, roc_node_2$tpr, roc_node_3$tpr)
thresholds <- list(
  roc_node_1$thresholds, roc_node_2$thresholds, roc_node_3$thresholds)
negative_count <- c(
  roc_node_1$negative_count, roc_node_2$negative_count, roc_node_3$negative_count)
total_count <- c(
  roc_node_1$total_count, roc_node_2$total_count, roc_node_3$total_count)

# Compute the global ROC curve for the model
roc_aggregated <- roc_curve(fpr, tpr, thresholds, negative_count, total_count)

# Calculate the AUC
roc_auc <- trapz(roc_aggregated$fpr, roc_aggregated$tpr)

sprintf("ROC AUC aggregated from each node's results: %f", roc_auc)

# Calculate the precision-recall
precision_recall_aggregated <- precision_recall_curve(
  fpr, tpr, thresholds, negative_count, total_count)

# Calculate the precision-recall AUC
precision_recall_auc <- -trapz(
  precision_recall_aggregated$recall, precision_recall_aggregated$pre)

sprintf(
  "Precision-Recall AUC aggregated from each node's results: %f",
  precision_recall_auc
)
```

## Validation

Using `ROCR` we can calculate the ROC and its AUC for the case of having all 
the data centrally available. The values between this and the aggregated ROC 
should match.

```{r validation}
roc_central_case <- get_roc(rbind(node_1, node_2, node_3))

# Validate the ROC AUC
sprintf(
  "ROC AUC using ROCR with all the data centrally available: %f",
  roc_central_case$auc@y.values[[1]]
)

# Validate the precision-recall AUC
precision_recall_auc <- trapz(
  roc_central_case$tpr,
  ifelse(is.nan(roc_central_case$prec), 1, roc_central_case$prec)
)
sprintf(
  "Precision-Recall AUC using ROCR with all the data centrally available: %f",
  precision_recall_auc
)
```

## Visualization

The ROC curve obtained can be visualized in the following way:

```{r visualization}
plot(roc_aggregated$fpr,
     roc_aggregated$tpr,
     main="ROC curve",
     xlab = "False Positive Rate",
     ylab = "True Positive Rate",
     cex=0.3,
     col="blue",
)
```

## Appendix

### Using pROC library

Another popular package to compute ROC curves is the `pROC`. Similarly to the 
example with the `ROCR` package, it's also possible to aggregate the results 
from ROC curves computed with the `pROC package`.

```{r appendix-proc}
library(pROC, warn.conflicts = FALSE)

get_proc <- function(dataset){
  glm.probs <- predict(glm.fit,
                       newdata = dataset,
                       type = "response")
  roc_obj <- roc(c(dataset$y), c(glm.probs))
  list(
    "fpr" = 1 - roc_obj$specificities,
    "tpr" = roc_obj$sensitivities,
    "thresholds" = roc_obj$thresholds,
    "negative_count"= sum(dataset$y == 0),
    "total_count" = nrow(dataset),
    "auc" = roc_obj$auc
  )
}

roc_obj_node_1 <- get_proc(node_1)
roc_obj_node_2 <- get_proc(node_2)
roc_obj_node_3 <- get_proc(node_3)

# Preparing the input
fpr <- list(roc_obj_node_1$fpr, roc_obj_node_2$fpr, roc_obj_node_3$fpr)
tpr <- list(roc_obj_node_1$tpr, roc_obj_node_2$tpr, roc_obj_node_3$tpr)
thresholds <- list(
  roc_obj_node_1$thresholds, roc_obj_node_2$thresholds, roc_obj_node_3$thresholds)
negative_count <- c(
  roc_obj_node_1$negative_count, roc_obj_node_2$negative_count, roc_obj_node_3$negative_count)
total_count <- c(
  roc_obj_node_1$total_count, roc_obj_node_2$total_count, roc_obj_node_3$total_count)

# Compute the global ROC curve for the model
roc_aggregated <- roc_curve(fpr, tpr, thresholds, negative_count, total_count)

# Calculate the AUC
roc_auc <- trapz(roc_aggregated$fpr, roc_aggregated$tpr)

sprintf("ROC AUC aggregated from each node's results: %f", roc_auc)

# Validate the ROC AUC
roc_central_case <- get_proc(rbind(node_1, node_2, node_3))

sprintf(
  "ROC AUC using pROC with all the data centrally available: %f",
  roc_central_case$auc
)
```
