---
title: "BY-COVID - WP5 - Baseline Use Case: SARS-CoV-2 vaccine effectiveness assessment"
subtitle: "Imputation of missing values"
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 - imputation of missing values'
  authors     <- 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 including the flag_violation_val (database table in BY-COVID-WP5-BaselineUseCase-VE.duckdb)')
  output     <- list('3_imputation.html (report)','cohort_data including the flag_violation_val and flag_listwise_del, cohort_data_imputed and imputation_methods (database tables in BY-COVID-WP5-BaselineUseCase-VE.duckdb)')
}

```

```{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 3_imputation.QMD
    
info(logger, paste0("
========================================================================================
",
"3_imputation.QDM","
",
"========================================================================================
 "))

```

```{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(
  {
    ## Description: get database connection
    con = dbConnect(duckdb::duckdb(), dbdir=auxiliary_database_path, read_only=FALSE)
    dbExecute(con, "CREATE OR REPLACE VIEW cohort_view AS SELECT * FROM cohort_data WHERE flag_violating_val==FALSE")
    df_not_viol <- dbGetQuery(con,"CREATE TEMPORARY TABLE Temp AS SELECT * FROM cohort_view;
    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)
  })

```

## Missing values {.tabset}

::: panel-tabset
##### Missing values map

```{r, missing values map, message=FALSE, warning=FALSE, error=FALSE}

info(logger_simple, "CREATING MISSING VALUES MAP...")

tryCatch(
  {
    explanatory <- c("age_nm","sex_cd","residence_area_cd","essential_worker_bl",
        "institutionalized_bl","foreign_bl", "diabetes_bl","obesity_bl","heart_failure_bl","copd_bl",
        "solid_tumor_without_metastasis_bl","chronic_liver_disease_bl","chronic_kidney_disease_bl",
        "sickle_cell_disease_bl","hypertension_bl","blood_cancer_bl","transplanted_bl", 
        "hiv_infection_bl","primary_immunodeficiency_bl","immunosuppression_bl","pregnancy_bl")
    pl_missing_map <- df_not_viol %>% missing_plot(explanatory = explanatory)
    ply <- ggplotly(pl_missing_map)
    for (i in 1:length(explanatory)) {
      ply$x$data[[1]]$text[i,] <- paste0("Id: ", ply$x$data[[1]][["x"]],"<br />","Variable: ", explanatory[length(explanatory)+1-i])
    }
    ply
  },
  error=function(cond) {
    ## Log info
    warn(logger, paste0("MY ERROR: 
                        ", cond))
    return(stop(cond))
  })

```

##### Missing data matrix

```{r, missing data matrix, message=FALSE, warning=FALSE, error=FALSE}

info(logger_simple, "CREATING MISSING DATA MATRIX...")

tryCatch(
  {
    explanatory <- c("age_nm","sex_cd","residence_area_cd","pregnancy_bl","essential_worker_bl",
        "institutionalized_bl","foreign_bl","socecon_lvl_cd") 
    df_not_viol %>% missing_pairs(explanatory = explanatory,position = "fill")
  },
  error=function(cond) {
    ## Log info
    warn(logger, paste0("MY ERROR: 
                        ", cond))
    return(stop(cond))
  })

```
:::

## Handling missing data {.tabset}

::: panel-tabset
##### Methods for handling missing data

```{r, variable status 1, results='hide', message=FALSE, warning=FALSE, error=FALSE}

## Description: create a table with variable characteristics

info(logger_simple, "CREATING A TABLE WITH VARIABLE CHARACTERISTICS...")

tryCatch(
  {
    ## Description: core variables (that require imputation of missing values)
    core_v <- 
      c("age_nm","sex_cd","residence_area_cd","pregnancy_bl","essential_worker_bl",
        "institutionalized_bl","foreign_bl","socecon_lvl_cd")
    
    ## Description: variables required to compute the variables 'comorbidities_bl' and 'immunestatus_bl'
    comorb_v <- c(
        ## comorbidities_bl:
        "diabetes_bl","obesity_bl","heart_failure_bl","copd_bl",
        "solid_tumor_without_metastasis_bl","chronic_liver_disease_bl","chronic_kidney_disease_bl",
        "sickle_cell_disease_bl","hypertension_bl")
    imm_v <- c(
        ## immunestatus_bl:
        "blood_cancer_bl","transplanted_bl",
        "hiv_infection_bl","primary_immunodeficiency_bl","immunosuppression_bl")
    
    ## Description: create missing data profile
    dt_prof <- profile_missing(df_not_viol)
    
    ## Description: create table to collect characteristics/checks for each variable
    df_var_imputation_method <- data.frame(Variable_name=dt_prof$feature)
    # Number and percentage missing for each variable
    df_var_imputation_method <- merge(x=df_var_imputation_method,y=dt_prof, 
          by.x=c("Variable_name"), 
          by.y=c("feature"),
          all.x=TRUE) 
    df_var_imputation_method <- df_var_imputation_method %>%
      mutate(Missing_values = ifelse(num_missing==0,FALSE,TRUE)) %>%
      mutate(Perc_miss_lt = ifelse(pct_missing>0.025,FALSE,TRUE)) %>%
      mutate(Perc_miss_lt5 = ifelse(pct_missing>0.05,FALSE,TRUE)) %>%
      mutate(Perc_miss_lt15 = ifelse(pct_missing>0.15,FALSE,TRUE)) # %>%
      # mutate(Perc_miss_lt50 = ifelse(pct_missing>0.50,FALSE,TRUE))
    df_var_imputation_method <- df_var_imputation_method %>% 
      mutate(Core = ifelse(Variable_name %in% core_v,TRUE,FALSE),
             Comorbidity = ifelse(Variable_name %in% comorb_v,TRUE,FALSE),
             Immunestatus = ifelse(Variable_name %in% imm_v,TRUE,FALSE))
    
    ## Description: check whether at least 3 'Comorbidity' variables less than 15% of values missing
    df_var_imputation_method$comorb_incl <- sum(df_var_imputation_method[which(df_var_imputation_method$Variable_name %in% comorb_v),"Perc_miss_lt15"])>=3
    ## Description: check whether at least 3 'Immune status' variables less than 15% of values missing
    df_var_imputation_method$imm_incl <- sum(df_var_imputation_method[which(df_var_imputation_method$Variable_name %in% imm_v),"Perc_miss_lt15"])>=3
  },
  error=function(cond) {
    ## Log info
    warn(logger, paste0("MY ERROR: 
                        ", cond))
    return(stop(cond))
  })
    
```

```{r, variable status 2, results='hide', message=FALSE, warning=FALSE, error=FALSE}

## Description: create a table with variable characteristics (add imputation method)

info(logger_simple, "CREATING A TABLE WITH VARIABLE CHARACTERISTICS (ADD IMPUTATION METHOD)...")

tryCatch(
  {
    ## Description: test MCAR assumption (exclude variables which have more than 50% of values missing)
    v_excl <- c()
    for(i in c(core_v, comorb_v, imm_v)) {
      if(length(unique(df_not_viol[[i]]))==1 |
         (df_var_imputation_method %>% filter(Variable_name==i) %>% select(Perc_miss_lt15))[[1]]==FALSE) {
        v_excl <- c(v_excl,i)
      }
    }
    mcar_res <- naniar::mcar_test(df_not_viol %>% select(all_of(setdiff(c(core_v,comorb_v, imm_v),v_excl)))) 
    # If p-value<0.05 reject the null hypothesis that the data are missing completely at random
    cond_mcar <- mcar_res$p.value>=0.05
    # cond_mcar <- FALSE
    df_var_imputation_method$MCAR <- cond_mcar
    
    ## Description: create methods table
    df_var_imputation_method <- df_var_imputation_method %>%
      mutate(imputation_method = case_when(
        Missing_values==FALSE ~ "No missing values",
        Missing_values==TRUE & Comorbidity==TRUE & comorb_incl==TRUE ~ "Allow missing values and don't impute missing values (comorbidity)",
        Missing_values==TRUE & Comorbidity==TRUE & comorb_incl==FALSE ~ "Exclude comorbidities_bl as matching variable (comorbidity)",
        Missing_values==TRUE & Immunestatus==TRUE & imm_incl==TRUE ~ "Allow missing values and don't impute missing values (immune status)",
        Missing_values==TRUE & Immunestatus==TRUE & imm_incl==FALSE ~ "Exclude immunestatus_bl as matching variable (immune status)",
        Missing_values==TRUE & Comorbidity==FALSE & Immunestatus==FALSE & Core==FALSE ~ "Allow missing values and don't impute missing values (no core variable)",
        Missing_values==TRUE & Comorbidity==FALSE & Immunestatus==FALSE & Core==TRUE & Perc_miss_lt15==FALSE ~ "Exclude core variable as matching variable (more than 15% missing values)",
        Missing_values==TRUE & Comorbidity==FALSE & Immunestatus==FALSE & Core==TRUE & Perc_miss_lt15==TRUE & cond_mcar==TRUE & Perc_miss_lt==TRUE ~ "Listwise deletion where core variable has missing values (MCAR reasonable)",
        Missing_values==TRUE & Comorbidity==FALSE & Immunestatus==FALSE & Core==TRUE & Perc_miss_lt15==TRUE & cond_mcar==TRUE & Perc_miss_lt==FALSE ~ "Imputation of missing values (MCAR reasonable)",
        Missing_values==TRUE & Comorbidity==FALSE & Immunestatus==FALSE & Core==TRUE & Perc_miss_lt15==TRUE & cond_mcar==FALSE & Perc_miss_lt5==TRUE ~ "Imputation of missing values (not MCAR)",
        Missing_values==TRUE & Comorbidity==FALSE & Immunestatus==FALSE & Core==TRUE & Perc_miss_lt15==TRUE & cond_mcar==FALSE & Perc_miss_lt5==FALSE ~ "Imputation of missing values Imputation of missing values (not MCAR, limit causal interpretation)"
      )) 
  },
  error=function(cond) {
    ## Log info
    warn(logger, paste0("MY ERROR: 
                        ", cond))
    return(stop(cond))
  })

```

```{r, save variable status, results='hide', message=FALSE, warning=FALSE, error=FALSE}

## Description: save table with variable characteristics

info(logger_simple, "SAVE TABLE WITH VARIABLE CHARACTERISTICS...")

tryCatch(
    {
      ## Description: Save table with variable characteristics as duckdb table
      con = dbConnect(duckdb::duckdb(), dbdir=auxiliary_database_path, read_only=FALSE)
      dbWriteTable(con, "imputation_methods",df_var_imputation_method,overwrite = TRUE, append=FALSE)
    },
    error=function(cond) {
      ## Log info
      warn(logger, paste0("MY ERROR: 
                          ", cond))
      return(stop(cond))
    },
    finally={
      dbDisconnect(con, shutdown=TRUE)
    })

```

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

## Description: listwise deletion

info(logger_simple, "LISTWISE DELETION...")

tryCatch(
  {
    ## Description: select variables that require listwise deletion
    v_ld <- as.character(df_var_imputation_method[which(df_var_imputation_method$imputation_method %in% c("Listwise deletion where core variable has missing values (MCAR reasonable)")),"Variable_name"])
  },
  error=function(cond) {
    ## Log info
    warn(logger, paste0("MY ERROR: 
                        ", cond))
    return(stop(cond))
  })

if(length(v_ld)>0) {
  tryCatch(
    {
      ## Description: variables that require listwise deletion
      tmp_ld_rn <- data.frame(person_id=df_not_viol[which(is.na(df_not_viol[,v_ld])),"person_id"])
      con = dbConnect(duckdb::duckdb(), dbdir=auxiliary_database_path, read_only=FALSE)
      duckdb_register(con, name='table_listwise_deletion', df=tmp_ld_rn, overwrite = TRUE, experimental = FALSE)
      qry <- 
        "update cohort_data set
          flag_listwise_del = TRUE WHERE cohort_data.person_id IN (select person_id from table_listwise_deletion);"
      dbExecute(con, qry)
      duckdb_unregister(con, name='table_listwise_deletion')
      rm(tmp_ld_rn)
    },
    error=function(cond) {
      ## Log info
      warn(logger, paste0("MY ERROR: 
                          ", cond))
      return(stop(cond))
    },
    finally={
      dbDisconnect(con, shutdown=TRUE)
    })
}

```

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

## Description: imputation of missing values

info(logger_simple, "IMPUTING MISSING VALUES...")

tryCatch(
  {
    ## Description: select variables that require imputation
    v_imp <- as.character(df_var_imputation_method[which(df_var_imputation_method$imputation_method %in% c("Imputation of missing values (MCAR reasonable)","Imputation of missing values (not MCAR)","Imputation of missing values (not MCAR, limit causal interpretation)")),"Variable_name"])
    
    ## Description: function for variables that require imputation
    method_imputation <- function(df_imp, v) {
      tmp <- df_imp %>% 
        select(-person_id) %>% 
        mutate_at(if("residence_area_cd" %in% names(.)) "residence_area_cd" else integer(0), factor) %>%
        mutate_at(if("sex_cd" %in% names(.)) "sex_cd" else integer(0), factor)

      # Imputation age: if sex_cd!="2" or age_nm<12 or age_nm>50 -> impute pregnancy_bl to false
      if("pregnancy_bl" %in% v) {
        tmp <- tmp %>% mutate(pregnancy_bl=ifelse(is.na(pregnancy_bl) & (sex_cd!="2" | age_nm<12 | age_nm>50),FALSE,pregnancy_bl))
      }
      if(length(v)>0) {
        ## If variables selected for imputation
        # -> Create prediction matrix
        pred <- make.predictorMatrix(tmp)
        # -> Adjust prediction matrix
        pred[,] <- 0
        for(j in v_imp) {
          pred[j,setdiff(c(core_v,comorb_v, imm_v),c(v_excl,j))] <- 1
        }
        # -> Imputation methods
        meth <- make.method(tmp)
        meth[] <- ""
        meth[v_imp] <- "pmm"
        # -> Run imputation
        imp <- mice(
          tmp,
          method = meth, 
          predictorMatrix = pred,
          m = 1, 
          maxit = 20
          )
        data_imp <- complete(imp) 
        data_imp <- cbind(person_id = df_imp$person_id, data_imp)
        # -> Select patient_ids of patients with imputed values
        out <- lapply(v, FUN =  function(variable) {
          out <- data.frame(person_id=df_imp[which(is.na(df_imp[,variable])),"person_id"])
          print(out)
          return(out)
        }) %>% bind_rows()
        p_id <- out %>% distinct()
        # -> Select records with imputed data and insert into duckdb
        data_imp <- data_imp[which(data_imp$person_id %in% p_id$person_id),]
      } else {
        ## If no variables selected for imputation
        data_imp <- df_imp[0,]
        imp <- NULL
      }
      return(list(data_imp,imp))
    }
    # Use the function to obtain an imputed dataset
    tmp <- method_imputation(df_not_viol %>% select(all_of(setdiff(c("person_id",core_v,comorb_v, imm_v),v_excl))),v_imp)
    data_imp <- tmp[[1]] %>% mutate(sex_cd=as.character(sex_cd), residence_area_cd=as.character(residence_area_cd))
  },
  error=function(cond) {
    ## Log info
    warn(logger, paste0("MY ERROR: 
                        ", cond))
    return(stop(cond))
  })

info(logger_simple, "CREATE DATABASE TABLE WITH IMPUTED DATA ONLY...")

tryCatch(
    {
      ## Create database table with imputed data only
      con = dbConnect(duckdb::duckdb(), dbdir=auxiliary_database_path, read_only=FALSE)
      dbWriteTable(con, "cohort_data_imputed",data_imp,overwrite = TRUE, append=FALSE)
      # print(dbGetQuery(conn = con, "SELECT * FROM cohort_data_imputed"))
    },
    error=function(cond) {
      ## Log info
      warn(logger, paste0("MY ERROR: 
                          ", cond))
      return(stop(cond))
    },
    finally={
      dbDisconnect(con, shutdown=TRUE)
    })

```

```{r, imputation process print, message=FALSE, warning=FALSE, error=FALSE}

## Description: create table with information on the imputation process

info(logger_simple, "CREATING A TABLE WITH INFORMATION ON THE IMPUTATION PROCESS...")

tryCatch(
  {
    ## Description: obtain number of imputed values
    df_var_imputation_method$Nr_imputed <- 0
    for(c in 1:length(v_imp)) {
      coln <- v_imp[c]
      df_var_imputation_method[which(df_var_imputation_method$Variable_name==coln),"Nr_imputed"] <- df_var_imputation_method[which(df_var_imputation_method$Variable_name==coln),"num_missing"]
    }
    df_var_imputation_method <- df_var_imputation_method %>% select(-c(num_missing,pct_missing))
    
    df_var_imputation_method <- df_var_imputation_method[, c("Variable_name","imputation_method","Nr_imputed","Missing_values","Comorbidity","comorb_incl","Immunestatus","imm_incl","Core","MCAR","Perc_miss_lt","Perc_miss_lt5","Perc_miss_lt15")]
    
    df_var_imputation_method <- df_var_imputation_method %>% arrange(factor(Variable_name,levels = colnames(df_not_viol)))
    
    ## Description: print table with information on methods used to deal with missing data
    formattable(df_var_imputation_method,
                align = c("l","c","c","c","c","c","c","c","c","c","c","c"),
                list(
                  `Method` = formatter("imputation_method",
                    style = ~style(display = "block", 
                                    padding = "0 4px", 
                                    `border-radius` = "4px", 
                                    `background-color` = case_when(`Method`=="No missing values" ~ "#E5FFCC", 
                                                                   `Method`=="Allow missing values and don't impute missing values (comorbidity)" ~ "#E0E0E0",
                                                                   `Method`=="Allow missing values and don't impute missing values (immune status)" ~ "#E0E0E0",
                                                                   `Method`=="Allow missing values and don't impute missing values (no core variable)" ~ "#E0E0E0",
                                                                   `Method`=="Imputation of missing values (MCAR reasonable)" ~"#F4A460",
                                                                   `Method`=="Imputation of missing values (not MCAR)" ~"#F4A460",
                                                                   `Method`=="Imputation of missing values (not MCAR, limit causal interpretation)" ~"#F4A460",
                                                                   `Method`=="Listwise deletion where core variable has missing values (MCAR reasonable)" ~ "#CCE5FF",
                                                                   `Method`=="Exclude comorbidities_bl as matching variable (comorbidity)" ~ "#cdccff",
                                                                   `Method`=="Exclude immunestatus_bl as matching variable (immune status)" ~ "#cdccff",
                                                                   `Method`=="Exclude core variable as matching variable (more than 15% missing values)" ~ "#cdccff"
                                                                   ))),
                  `Number of imputed values` = formatter("Nr_imputed",
                    style = ~style(display = "block", 
                                    padding = "0 4px", 
                                    `border-radius` = "4px", 
                                    `background-color` = case_when(`Number of imputed values`>0 ~"#F4A460")))
                  )
                ) %>% rename("Variable" = Variable_name,"Method" = imputation_method, "Number of imputed values" = Nr_imputed)
  },
  error=function(cond) {
    ## Log info
    warn(logger, paste0("MY ERROR: 
                        ", cond))
    return(stop(cond))
  })

```

##### Distribution imputed data

```{r, plot imputation, message=FALSE, warning=FALSE, error=FALSE}

## Description: plot the distribution of imputed data

info(logger_simple, "PLOT THE DISTRIBUTION OF IMPUTED DATA...")

tryCatch(
  {
    imp <- tmp[[2]]
    
    if(!is.null(imp)) {
      ## Create plots
      plot(imp)
      formula <- as.formula(paste("~", paste(v_imp, 
            collapse = "+", sep = ""), sep = ""))
      densityplot(imp,formula)
    }
  },
  error=function(cond) {
    ## Log info
    warn(logger, paste0("MY ERROR: 
                        ", cond))
    return(stop(cond))
  })

```
:::

## Listwise deletion

```{r, print listwise deletion, message=FALSE, warning=FALSE, error=FALSE}

## Description: select patient ids of records not used for further analysis because of missing values (listwise deletion)

info(logger_simple, "ADJUSTING FLAG VARIABLE (FLAG_LISTWISE_DEL)...
     ")

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

```
