---
title: "CohortConstructor benchmarking results"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{a11_benchmark}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---


```{r, include = FALSE}
NOT_CRAN <- identical(tolower(Sys.getenv("NOT_CRAN")), "true")

knitr::opts_chunk$set(
collapse = TRUE,
eval = TRUE, 
warning = FALSE, 
message = FALSE,
comment = "#>",
echo = FALSE,
eval = NOT_CRAN
)
```
  
```{r, include = FALSE}
library(CDMConnector)
if (Sys.getenv("EUNOMIA_DATA_FOLDER") == ""){
Sys.setenv("EUNOMIA_DATA_FOLDER" = file.path(tempdir(), "eunomia"))}
if (!dir.exists(Sys.getenv("EUNOMIA_DATA_FOLDER"))){ dir.create(Sys.getenv("EUNOMIA_DATA_FOLDER"))
downloadEunomiaData()  
}
```


```{r}
# Packages
library(visOmopResults)
library(omopgenerics)
library(ggplot2)
library(CohortCharacteristics)
library(stringr)
library(dplyr)
library(tidyr)
library(gt)
library(scales)
library(CohortConstructor)
library(gt)

niceOverlapLabels <- function(labels) {
  new_labels <- gsub("_", " ", gsub(" and.*|cc_", "", labels))
  return(
    tibble("Cohort name" = new_labels) |>
      mutate(
        "Cohort name" = str_to_sentence(gsub("_", " ", gsub("cc_|atlas_", "", new_labels))),
        "Cohort name" = case_when(
          grepl("Asthma", .data[["Cohort name"]]) ~ "Asthma without COPD",
          grepl("Covid", .data[["Cohort name"]]) ~ gsub("Covid|Covid", "COVID-19", `Cohort name`),
          grepl("eutropenia", .data[["Cohort name"]]) ~ "Acquired neutropenia or unspecified leukopenia",
          grepl("Hosp", .data[["Cohort name"]]) ~ "Inpatient hospitalisation",
          grepl("First", .data[["Cohort name"]]) ~ "First major depression",
          grepl("fluoro", .data[["Cohort name"]]) ~ "New fluoroquinolone users",
          grepl("Beta", .data[["Cohort name"]]) ~ "New users of beta blockers nested in essential hypertension",
          .default = .data[["Cohort name"]]
        ),
        "Cohort name" = if_else(
          grepl("COVID", .data[["Cohort name"]]),
          gsub(" female", ": female", gsub(" male", ": male", .data[["Cohort name"]])),
          .data[["Cohort name"]]
        ),
        "Cohort name" = if_else(
          grepl(" to ", .data[["Cohort name"]]),
          gsub("male ", "male, ", .data[["Cohort name"]]),
          .data[["Cohort name"]]
        )
      )
  )
}
```

# Introduction

Cohorts are a fundamental building block for studies that use the OMOP CDM, identifying people who satisfy one or more inclusion criteria for a duration of time based on their clinical records. Currently cohorts are typically built using [CIRCE](https://github.com/OHDSI/circe-be) which allows complex cohorts to be represented using JSON. This JSON is then converted to SQL for execution against a database containing data mapped to the OMOP CDM. CIRCE JSON can be created via the [ATLAS](https://github.com/OHDSI/Atlas) GUI or programmatically via the [Capr](https://github.com/OHDSI/Capr) R package. However, although a powerful tool for expressing and operationalising cohort definitions, the SQL generated can be cumbersome especially for complex cohort definitions, moreover cohorts are instantiated independently, leading to duplicated work.

The CohortConstructor package offers an alternative approach, emphasising cohort building in a pipeline format. It first creates base cohorts and then applies specific inclusion criteria. Unlike the "by definition" approach, where cohorts are built independently, CohortConstructor follows a "by domain/ table" approach, which minimises redundant queries to large OMOP tables. More details on this approach can be found in the [Introduction vignette](https://ohdsi.github.io/CohortConstructor/articles/a00_introduction.html).

To test the performance of the package there is a benchmarking function which uses nine phenotypes from the [OHDSI Phenotype library](https://github.com/OHDSI/PhenotypeLibrary) that cover a range of concept domains, entry and inclusion criteria, and cohort exit options. We replicated these cohorts using CohortConstructor to assess computational time and agreement between CIRCE and CohortConstructor.

```{r, echo=TRUE}
library(CDMConnector)
library(CodelistGenerator)
library(PatientProfiles)
library(CohortConstructor)
library(dplyr)

con <- DBI::dbConnect(duckdb::duckdb(), 
                      dbdir = eunomiaDir())
cdm <- cdmFromCon(con, cdmSchema = "main", writeSchema = "main", 
                  writePrefix = "my_study_")
```

Once we have created our cdm reference we can run the benchmark. Once run we'll have a set of results with the time taken to run the different tasks. For this example we will just run task of creating all the cohorts at once using CohortConstructor.

```{r, echo=TRUE}
benchmark_results <- benchmarkCohortConstructor(
  cdm,
  runCIRCE = FALSE,
  runCohortConstructorDefinition = FALSE,
  runCohortConstructorDomain = TRUE
)
benchmark_results |> 
  glimpse()
```

## Collaboration

If you are interested in running the code on your database, feel free to reach out to us for assistance, and we can also update the vignette with your results! :)

The benchmark script was executed against the following databases:

- **CPRD Gold**: A primary care database from the UK, capturing data mostly from Northern Ireland, Wales, and Scotland clinics. The benchmark utilized a 100,000-person sample from this dataset, which is managed using PostgreSQL.

- **CPRD Aurum**: Another UK primary care database, primarily covering clinics in England. This database is managed on SQL Server.

- **Coriva**: A sample of approximately 400,000 patients from the Estonia National Health Insurance database, managed on PostgreSQL.

- **OHDSI SQL Server**: A mock OMOP CDM dataset provided by OHDSI, hosted on SQL Server.

The table below presents the number of records in the OMOP tables used in the benchmark script for each of the participating databases.

```{r}
benchmarkData$omop |>
  visOmopResults::formatTable() |>
  tab_style(style = list(cell_fill(color = "#e1e1e1"), cell_text(weight = "bold")),
            locations = cells_column_labels()) |>
  tab_style(style = list(cell_text(weight = "bold")),
            locations = cells_body(columns = 1))
```



# Cohorts

We replicated the following cohorts from the OHDSI phenotype library: COVID-19 (ID 56), inpatient hospitalisation (23),  new users of beta blockers nested in essential hypertension (1049), transverse myelitis (63), major non cardiac surgery (1289), asthma without COPD (27), endometriosis procedure (722), new fluoroquinolone users (1043), acquired neutropenia or unspecified leukopenia (213). 

The COVID-19 cohort was used to evaluate the performance of common cohort stratifications. To compare the package with CIRCE, we created definitions in Atlas, stratified by age groups and sex, which are available in the [benchmark GitHub repository](https://github.com/oxford-pharmacoepi/BenchmarkCohortConstructor/tree/main/JSONCohorts) with the benchmark code.

## Cohort counts and overlap

The following table displays the number of records and subjects for each cohort across the participating databases:

```{r}
benchmarkData$details |>
  visOmopResults::formatTable(groupColumn = "cdm_name") |>
  tab_style(style = list(cell_fill(color = "#e1e1e1"), cell_text(weight = "bold")),
            locations = cells_column_labels()) |>
  tab_style(style = list(cell_text(weight = "bold")),
            locations = cells_body(columns = 1:2))
```

We also computed the overlap between patients in CIRCE and CohortConstructor cohorts, with results shown in the plot below:

```{r, fig.width=10, fig.height=7}
gg <- benchmarkData$comparison |>
  plotCohortOverlap(uniqueCombinations = FALSE, facet = "cdm_name") +
  scale_y_discrete(labels = niceOverlapLabels) +
  theme(
    legend.text = element_text(size = 10),
    strip.text = element_text(size = 14),
    axis.text.x = element_text(size = 12),
    axis.title.x = element_text(size = 14),
    axis.title.y = element_text(size = 14)
  ) +
  # facet_wrap("cdm_name") +
  scale_fill_discrete(labels = c("Both", "CIRCE", "CohortConstructor")) +
  scale_color_discrete(labels = c("Both", "CIRCE", "CohortConstructor")) +
  xlab("") 

gg$data <- gg$data |>
  dplyr::mutate(
    cohort_name_reference_cohort_name_comparator = visOmopResults::customiseText(gsub("cc_| -.*", "", cohort_name_reference_cohort_name_comparator)),
    cohort_name_reference_cohort_name_comparator = case_when(
      grepl("Asthma", .data[["cohort_name_reference_cohort_name_comparator"]]) ~ "Asthma without COPD",
      grepl("Covid", .data[["cohort_name_reference_cohort_name_comparator"]]) ~ gsub("Covid|Covid", "COVID-19", cohort_name_reference_cohort_name_comparator),
      grepl("eutropenia", .data[["cohort_name_reference_cohort_name_comparator"]]) ~ "Acquired neutropenia or unspecified leukopenia",
      grepl("Hosp", .data[["cohort_name_reference_cohort_name_comparator"]]) ~ "Inpatient hospitalisation",
      grepl("First", .data[["cohort_name_reference_cohort_name_comparator"]]) ~ "First major depression",
      grepl("fluoro", .data[["cohort_name_reference_cohort_name_comparator"]]) ~ "New fluoroquinolone users",
      grepl("Beta", .data[["cohort_name_reference_cohort_name_comparator"]]) ~ "New users of beta blockers nested in essential hypertension",
      .default = .data[["cohort_name_reference_cohort_name_comparator"]]
    ),
    cohort_name_reference_cohort_name_comparator = if_else(
      grepl("COVID", .data[["cohort_name_reference_cohort_name_comparator"]]),
      gsub(" female", ": female", gsub(" male", ": male", .data[["cohort_name_reference_cohort_name_comparator"]])),
      .data[["cohort_name_reference_cohort_name_comparator"]]
    ),
    cohort_name_reference_cohort_name_comparator = if_else(
      grepl(" to ", .data[["cohort_name_reference_cohort_name_comparator"]]),
      gsub("male ", "male, ", .data[["cohort_name_reference_cohort_name_comparator"]]),
      .data[["cohort_name_reference_cohort_name_comparator"]]
    )
  )
gg
```

# Performance

To evaluate CohortConstructor performance we generated each of the CIRCE cohorts using functionalities provided by both CodelistGenerator and CohortConstructor, and measured the computational time taken.

Two different approaches with CohortConstructor were tested: 

- *By definition*: we created each of the cohorts seprately.

- *By domain*: All nine targeted cohorts were created together in a set, following the by domain approach described in the [Introduction vignette](https://ohdsi.github.io/CohortConstructor/articles/a00_introduction.html). Briefly, this approach involves creating all base cohorts at once, requiring only one call to each involved OMOP table.


## By definition

The following plot shows the times taken to create each cohort using CIRCE and CohortConstructor when each cohorts were created separately.

```{r}
## TABLE with same results as the plot below.

# header_prefix <- "[header]Time by database (minutes)\n[header_level]"
# benchmarkData$time |>
#   distinct() |>
#   filter(!grepl("male|set", msg)) |>
#   mutate(
#     time = niceNum((as.numeric(toc) - as.numeric(tic))/60, 2),
#     Tool = if_else(grepl("cc", msg), "CohortConstructor", "CIRCE"),
#     "Cohort name" = str_to_sentence(gsub("_", " ", gsub("cc_|atlas_", "", msg)))
#   ) |>
#   select(-c("tic", "toc", "msg", "callback_msg")) |>
#   pivot_wider(names_from = "cdm_name", values_from = "time", names_prefix = header_prefix) |>
#   select(c("Cohort name", "Tool", paste0(header_prefix, data$time$cdm_name |> unique()))) |>
#   mutate(
#     "Cohort name" = case_when(
#       grepl("Asthma", .data[["Cohort name"]]) ~ "Asthma without COPD",
#       grepl("Covid", .data[["Cohort name"]]) ~ "COVID-19",
#       grepl("eutropenia", .data[["Cohort name"]]) ~ "Acquired neutropenia or unspecified leukopenia",
#       grepl("Hosp", .data[["Cohort name"]]) ~ "Inpatient hospitalisation",
#       grepl("First", .data[["Cohort name"]]) ~ "First major depression",
#       grepl("fluoro", .data[["Cohort name"]]) ~ "New fluoroquinolone users",
#       grepl("Beta", .data[["Cohort name"]]) ~ "New users of beta blockers nested in essential hypertension",
#       .default = .data[["Cohort name"]]
#     )
#   ) |>
#   arrange(`Cohort name`) |>
#   gtTable(colsToMergeRows = "all_columns") |>
#   tab_style(style = list(cell_fill(color = "#e1e1e1"), cell_text(weight = "bold")), 
#             locations = cells_column_labels()) |>
#   tab_style(style = list(cell_text(weight = "bold")), 
#             locations = cells_body(columns = 1:2)) 
```

```{r, fig.width=10, fig.height=7}

benchmarkData$time_definition |>
  ggplot(aes(y = `Cohort name`, x = time, colour = Tool, fill = Tool)) +
  geom_col(position = "dodge", width = 0.6) +
  xlab("Time (minutes)") +
  scale_y_discrete(labels = label_wrap(20)) +
  theme(
    legend.title = element_blank(),
    legend.position = "bottom",
    axis.text.x = element_text(size = 12),
    legend.text = element_text(size = 12),
    strip.text = element_text(size = 14),
    axis.text.y = element_text(size = 12),  
    axis.title.x = element_text(size = 14),  
    axis.title.y = element_text(size = 14)   
  ) +
  facet_wrap(vars(cdm_name), nrow = 1, scales = "free_x")
```

## By domain

The table below depicts the total time it took to create the nine cohorts when using the *by domain* approach for CohortConstructor.

```{r}
header_prefix <- "[header]Time by tool (minutes)\n[header_level]"
benchmarkData$time_domain |>
  formatTable() |> 
  tab_style(style = list(cell_fill(color = "#e1e1e1"), cell_text(weight = "bold")), 
            locations = cells_column_labels()) |>
  tab_style(style = list(cell_text(weight = "bold")), 
            locations = cells_body(columns = 1)) 
```

## Cohort stratification

Cohorts are often stratified in studies. With Atlas cohort definitions, each stratum requires a new CIRCE JSON to be instantiated, while CohortConstructor allows stratifications to be generated from an overall cohort. The following table shows the time taken to create age and sex stratifications for the COVID-19 cohort with both CIRCE and CohortConstructor.

```{r}
benchmarkData$time_strata |>
  formatTable() |>
  tab_style(style = list(cell_fill(color = "#e1e1e1"), cell_text(weight = "bold")), 
            locations = cells_column_labels()) |>
  tab_style(style = list(cell_text(weight = "bold")), 
            locations = cells_body(columns = 1)) 
```
