--- 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 ) ```