---
title: "BY-COVID - WP5 - Baseline Use Case: SARS-CoV-2 vaccine effectiveness assessment"
subtitle: "Validation"
format: 
  html:
    self-contained: true
code-fold: false
echo: false
editor: visual
always_allow_html: yes
title-block-banner: "#27445C"
---

```{r, metadata, results='hide', message=FALSE, warning=FALSE, error=FALSE}

################
### Metadata ###
################

# DATE LAST MODIFIED:
# 19/09/2023

# METADATA: 
if(FALSE) {
  title      <- 'BY-COVID WP5.2 Baseline Use Case: SARS-CoV-2 vaccine effectiveness - analytical pipeline - validation'
  author     <- list('Marjan Meurisse','Javier González-Galindo','Francisco Estupiñán-Romero','Santiago Royo-Sierra','Nina Van Goethem','Enrique Bernal-Delgado')
  version    <- '1.0.2'
  maintainer <- 'Marjan Meurisse'
  email      <- 'Marjan.Meurisse@sciensano.be'
  input      <- list('cohort_data (database table in BY-COVID-WP5-BaselineUseCase-VE.duckdb)')
  output     <- list('2_validation.html (report)','Updated cohort_data table including the flag_violation_val')
}

```

```{css}

.box {
  padding: 1em;
  color: black;
  background: transparant;
  border: 1px solid lightgrey;
}

```

```{r, general settings, results='hide', message=FALSE, warning=FALSE, error=FALSE}

## Description: load packages
source("./0_global.R")

```

```{r, logging, results='hide', message=FALSE, warning=FALSE, error=FALSE}

## Description: Log of 2_validation.QMD
    
info(logger, paste0("
========================================================================================
",
"2_validation.QMD","
",
"========================================================================================
 "))

```

```{r, load data, results='hide', message=FALSE, warning=FALSE, error=FALSE}

## Description: load the data into memory

info(logger_simple, "LOADING THE DATA INTO MEMORY...")

tryCatch(
  {
    con = dbConnect(duckdb::duckdb(), dbdir=auxiliary_database_path, read_only=FALSE)
    df <- dbGetQuery(con,
    "CREATE TEMPORARY TABLE Temp AS SELECT * FROM cohort_data;
    ALTER TABLE Temp DROP COLUMN flag_violating_val;
    ALTER TABLE Temp DROP COLUMN flag_listwise_del;
    SELECT * FROM Temp")
  },
  error=function(cond) {
    ## Log info
    warn(logger, paste0("MY ERROR: 
                        ", cond))
    return(stop(cond))
  },
  finally={
    dbDisconnect(con, shutdown=TRUE)
  })

```

## Compliance with the Common Data Model specification {.tabset}

We check whether the imported dataset complies with the data model specification (<https://doi.org/10.5281/zenodo.6913045>).

To comply with the data model, the dataset must pass a number of validation rules. The data are tested against this set of validation rules and results from this validation process are summarized.

::: panel-tabset
```{r, validation rules, message=FALSE, warning=FALSE, error=FALSE}

## Description: specification of validation rules and confrontation of the data with these validation rules

info(logger_simple, "CHECKING DATA AGAINST A SET OF VALIDATION RULES...")

validation <- function(data) {
  
  ## specify validation rules
  rules <- validator(
    is.na(age_nm) | age_nm >= 5 & age_nm <=115,
    is.na(sex_cd) | sex_cd %in% c("0","1","2","9"),
    is.na(dose_1_brand_cd) | dose_1_brand_cd %in% c("BP","MD","JJ","AZ","NV"),
    is.na(dose_2_brand_cd) | dose_2_brand_cd %in% c("BP","MD","JJ","AZ","NV"),
    is.na(doses_nm) | doses_nm >= 0 & doses_nm <= 10,
    fully_vaccinated_bl==FALSE |
      fully_vaccinated_bl==TRUE & !is.na(vaccination_schedule_cd),
    is.na(test_type_cd) | test_type_cd %in% c("PCR","AG","other"),
    is.na(variant_cd) | variant_cd %in% c("alpha","beta", "gamma", "delta", "omicron","epsilon", "zeta", "eta", "theta", "iota", "kappa", "lambda", "mu"),
    is.na(pregnancy_bl) | pregnancy_bl==FALSE | (pregnancy_bl==TRUE & sex_cd=="2" & age_nm>=12 & age_nm<=55),
    is.na(essential_worker_bl) | essential_worker_bl==FALSE | (essential_worker_bl==TRUE & age_nm>=16 & age_nm<=70),
    
    (is.na(dose_1_dt) & is.na(dose_2_dt)) | is.na(dose_2_dt) |
      !is.na(dose_1_dt) & !is.na(dose_2_dt) & 
      (dose_1_dt < dose_2_dt),
    (is.na(dose_2_dt) & is.na(dose_3_dt)) | is.na(dose_3_dt) |
      !is.na(dose_2_dt) & !is.na(dose_3_dt) & 
      (dose_2_dt < dose_3_dt),
    is.na(previous_infection_dt) | is.na(confirmed_case_dt) |
      !is.na(previous_infection_dt) & !is.na(confirmed_case_dt) & 
      (previous_infection_dt < confirmed_case_dt),
    is.na(confirmed_case_dt) | is.na(exitus_dt) |
      !is.na(confirmed_case_dt) & !is.na(exitus_dt) & 
      (confirmed_case_dt <= exitus_dt),
    is.na(previous_infection_dt) | is.na(exitus_dt) |
      !is.na(previous_infection_dt) & !is.na(exitus_dt) & 
      (previous_infection_dt <= exitus_dt),
    is.na(fully_vaccinated_dt) | is.na(exitus_dt) |
      !is.na(fully_vaccinated_dt) & !is.na(exitus_dt) & 
      fully_vaccinated_dt <= exitus_dt,
    
    (!is.na(dose_1_dt) & !is.na(dose_2_dt) & !is.na(dose_3_dt) & doses_nm>=3) | 
        (!is.na(dose_1_dt) & !is.na(dose_2_dt) & is.na(dose_3_dt) & doses_nm==2) |
        (!is.na(dose_1_dt) & is.na(dose_2_dt) & is.na(dose_3_dt) & doses_nm==1) | 
        (is.na(dose_1_dt) & is.na(dose_2_dt) & is.na(dose_3_dt) & doses_nm==0),
    is.na(dose_1_dt) | (!is.na(dose_1_dt) & !is.na(dose_1_brand_cd)),
    is.na(dose_2_dt) | (!is.na(dose_2_dt) & !is.na(dose_2_brand_cd) &
                        !is.na(dose_1_dt) & !is.na(dose_1_brand_cd)),
    is.na(dose_3_dt) | (!is.na(dose_3_dt) & !is.na(dose_3_brand_cd) &
                        !is.na(dose_2_dt) & !is.na(dose_2_brand_cd) &
                        !is.na(dose_1_dt) & !is.na(dose_1_brand_cd)),
    (dose_1_brand_cd=="JJ" & !is.na(dose_1_dt) & !is.na(fully_vaccinated_dt) & fully_vaccinated_bl==TRUE) | 
      (dose_1_brand_cd!="JJ" & !is.na(dose_2_dt) & !is.na(fully_vaccinated_dt) & fully_vaccinated_bl==TRUE) | 
      (is.na(dose_1_brand_cd) & is.na(fully_vaccinated_dt) & fully_vaccinated_bl==FALSE) |
      (dose_1_brand_cd!="JJ" & is.na(dose_2_dt) & is.na(fully_vaccinated_dt) & fully_vaccinated_bl==FALSE)
  )
  
  ## confront the data with the validation rules
  out <- confront(data,rules)
  return(out)
}

tryCatch(
  {
    out <- validation(df)
    summary_validation_tmp <- summary(out)
  },
  error=function(cond) {
    ## Log info
    warn(logger, paste0("MY ERROR: 
                        ", cond))
    return(stop(cond))
  })

```

##### Validation table

```{r, validation table, message=FALSE, warning=FALSE, error=FALSE}

## Description: create a table with a summary of the validation results

info(logger_simple, "CREATE TABLE WITH SUMMARY OF VALIDATION RESULTS...")

tryCatch(
  {
    summary_validation <- summary_validation_tmp
    summary_validation$perc_fails <- paste0(round(100*summary_validation$fails/summary_validation$items,2),'%')
    summary_validation$perc_missing <- paste0(round(100*summary_validation$nNA/summary_validation$items,2),'%')
    summary_validation$perc_missing_num <- 100*summary_validation$nNA/summary_validation$items
    summary_validation <- summary_validation[,c("expression","name", "items","passes","fails","perc_fails","nNA","perc_missing", "error","warning","perc_missing_num")]
    colnames(summary_validation) <- c("Validation rule","Name rule","Items","Passes","Fails","Percentage of fails","Number of NAs","Percentage of NAs","Error","Warning","perc_missing_num")
    formattable(summary_validation,
                align = c("l","c","c","c","c","c","c","c","c","c","c"),
                list(
                  perc_missing_num = FALSE,
                  Error = formatter("Error",
                    style = ~style(display = "block", 
                                    padding = "0 4px", 
                                    `border-radius` = "4px", 
                                    `background-color` = case_when(Error==TRUE ~"#F4A460"))),
                  Warning = formatter("Warning",
                    style = ~style(display = "block", 
                                    padding = "0 4px", 
                                    `border-radius` = "4px", 
                                    `background-color` = case_when(Warning==TRUE ~"#F4A460"))),
                  `Percentage of NAs` = formatter("Percentage of NAs",
                    style = ~style(display = "block", 
                                    padding = "0 4px", 
                                    `border-radius` = "4px", 
                                    `background-color` = 
                                    case_when(perc_missing_num > 2.5 ~ "#F4A460"))),
                  Fails = formatter("Fails",
                    style = ~style(display = "block", 
                                    padding = "0 4px", 
                                    `border-radius` = "4px", 
                                    `background-color` = case_when(Fails > 0 ~ "#F4A460")))
    ))
  },
  error=function(cond) {
    ## Log info
    warn(logger, paste0("MY ERROR: 
                        ", cond))
    return(stop(cond))
  })

```

##### Validation plot

The vertical bars in the validation plot indicate the percentage of records 'Passing', 'Failing' and 'Missing'

```{r, validation plot 1, results='hide', fig.show='hide', message=FALSE, warning=FALSE, error=FALSE}

## Description: create interactive validation plot

info(logger_simple, "CREATE VALIDATION PLOT...")

tryCatch(
  {
    ## Create plot
    pl <- summary_validation %>%
      mutate(perc_fails=Fails/Items) %>%
      mutate(perc_passes=Passes/Items) %>%
      mutate(perc_missing=`Number of NAs`/Items) %>%
      select(c(`Name rule`,`Validation rule`,perc_fails,perc_passes,perc_missing)) %>%
      gather(name, value, perc_fails:perc_missing) %>%
      mutate(name=case_when(
        name=="perc_fails" ~ "Fails",
        name=="perc_passes" ~ "Passes",
        name=="perc_missing" ~ "Missing"
      )) %>%
      mutate(name=factor(name,levels=c("Fails","Passes","Missing"))) %>%
      ggplot(aes(x=`Name rule`,y=value,fill=name,order=name)) +
        geom_col(position = position_stack(reverse = TRUE)) +
        scale_y_continuous(labels=scales::percent) +
        theme_minimal() +
        theme(axis.title.x =element_blank(),
              axis.title.y = element_blank(),
              axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
              legend.title = element_blank()) +
        scale_fill_manual(
          name ="",
          values = c("Fails" = "red",
                  "Passes" = "#4C9900",
                  "Missing" = "grey"),
          breaks=c("Fails","Passes","Missing")) 
    
    ## Convert to plotly
    ply <- ggplotly(pl)
    
    ## Adjust hover text
    for (i in 1:length(ply$x$data)) {
      ply$x$data[[i]]$text <- paste0(gsub("<.*$", "", ply$x$data[[i]]$text), "<br />", "Value: ", ply$x$data[[i]]$y, "<br />", ply$x$data[[i]]$name)
    }
  },
  error=function(cond) {
    ## Log info
    warn(logger, paste0("MY ERROR: 
                        ", cond))
    return(stop(cond))
  })

```

```{r, validation plot 2, message=FALSE, warning=FALSE, error=FALSE}

tryCatch(
  {
    ply
  },
  error=function(cond) {
    ## Log info
    warn(logger, paste0("MY ERROR: 
                        ", cond))
    return(stop(cond))
  })

```
:::

## Non-compliance with the Common Data Model specification

The set of validation rules are considered 'essential' not to be violated to be considered for the subsequent analysis. A logical variable *flag_violation_val* is created in the cohort_data table in the BY-COVID-WP5-BaselineUseCase-VE.duckdb database and set to TRUE when at least one of the validation rules in the pre-specified set is violated (otherwise this variable is set to FALSE).

```{r, violating essential rules, message=FALSE, warning=FALSE, error=FALSE}

## Description: select patient ids of records violating at least one of the validation rules

info(logger_simple, "SELECTING PATIENT IDS OF RECORDS VIOLATING AT LEAST ONE OF THE VALIDATION RULES...")

tryCatch(
  {
    tmp_viol <- data.frame(person_id=violating(df,out)[,"person_id"])
  },
  error=function(cond) {
    ## Log info
    warn(logger, paste0("MY ERROR: 
                        ", cond))
    return(stop(cond))
  })

```

```{r, save violated, message=FALSE, warning=FALSE, error=FALSE}

## Description: flag_violating_val to TRUE when at least one of the rules is violated

info(logger_simple, "SET FLAG_VIOLATING_VAL TO TRUE WHEN AT LEAST ONE OF THE RULES IS VIOLATED...")

tryCatch(
  {
    con = dbConnect(duckdb::duckdb(), dbdir=auxiliary_database_path, read_only=FALSE)
    duckdb_register(con, name='table_violating_ids', df=tmp_viol, overwrite = TRUE, experimental = FALSE)
    qry <- 
      "update cohort_data set
        flag_violating_val = TRUE WHERE cohort_data.person_id IN (select person_id from table_violating_ids);"
    dbExecute(con, qry)
    duckdb_unregister(con, name='table_violating_ids')
  },
  error=function(cond) {
    ## Log info
    warn(logger, paste0("MY ERROR: 
                        ", cond))
    return(stop(cond))
  },
  finally={
    dbDisconnect(con, shutdown=TRUE)
  })

```

```{r, print violating rules, message=FALSE, warning=FALSE, error=FALSE}

## Description: count patients with flag_violating_val TRUE versus FALSE

info(logger_simple, "COUNT PATIENTS WITH FLAG_VIOLATING_VAL TRUE VERSUS FALSE...
     ")

tryCatch(
  {
    con = dbConnect(duckdb::duckdb(), dbdir=auxiliary_database_path, read_only=FALSE) 
    count_viol_true <- dbGetQuery(conn = con, "SELECT COUNT(*) FROM cohort_data WHERE flag_violating_val==TRUE")
    count_viol_false <- dbGetQuery(conn = con, "SELECT COUNT(*) FROM cohort_data WHERE flag_violating_val==FALSE")
    count_viol <- data.frame(matrix(ncol = 2, nrow = 0))
    count_viol <- rbind(count_viol,c(count_viol_true,count_viol_false))
    colnames(count_viol) <- c("flag_violating_val==TRUE","flag_violating_val==FALSE")
  },
  error=function(cond) {
    ## Log info
    warn(logger, paste0("MY ERROR: 
                        ", cond))
    return(stop(cond))
  },
  finally={
    dbDisconnect(con, shutdown=TRUE)
  })

formattable(count_viol,
            align="c")

```
