---
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:
# 20/04/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.0'
  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, and cohort_data_imputed (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, load data, results='hide', message=FALSE, warning=FALSE, error=FALSE}

## Description: Load the data into memory

tryCatch(
  {
    con = dbConnect(duckdb::duckdb(), dbdir=auxilary_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")
  },
  finally={
    dbDisconnect(con, shutdown=TRUE)
  })

```

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

## Description: variables indicated as required in the data model
required_v <- 
  c("person_id","age_nm","sex_cd","residence_area_cd","country_cd","foreign_bl",
    "exitus_dt","exitus_bl",
    "essential_worker_bl","institutionalized_bl", 
    "dose_1_brand_cd", "dose_1_dt","dose_2_brand_cd","dose_2_dt",
    "fully_vaccinated_dt","fully_vaccinated_bl","vaccination_schedule_cd",
    "confirmed_case_dt","confirmed_case_bl","previous_infection_dt","previous_infection_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")

## 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")

## Description: variables required to compute the variables 'comorbidities_bl' and 'immunestatus_bl'
comorb_imm_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",
    ## 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))
df_var_imputation_method <- df_var_imputation_method %>%
  mutate(All_missing_values = ifelse(pct_missing==1,TRUE,FALSE))
df_var_imputation_method <- df_var_imputation_method %>% 
  mutate(Required = ifelse(Variable_name %in% required_v,TRUE,FALSE)) %>%
  mutate(Core = ifelse(Variable_name %in% core_v,TRUE,FALSE)) %>%
  mutate(Comorbidity = ifelse(Variable_name %in% comorb_imm_v,TRUE,FALSE)) 

```

## Imputation {.tabset}

::: panel-tabset
##### Imputation process

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

## Description: test MCAR assumption

v_equal <- c()
for(i in core_v) {
  if(length(unique(df_not_viol[[i]]))==1) {
    v_equal <- c(v_equal,i)
  }
}
ifelse(df_var_imputation_method[which(df_var_imputation_method$Variable_name=="pregnancy_bl"),"All_missing_values"],
       mcar_res <- naniar::mcar_test(df_not_viol %>% select(any_of(setdiff(colnames(df_not_viol),v_equal))) %>% select(all_of(setdiff(core_v,"pregnancy_bl")))),
       mcar_res <- naniar::mcar_test(df_not_viol %>% select(any_of(setdiff(colnames(df_not_viol),v_equal))) %>% select(all_of(core_v)))   
       )
# 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
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 & Required==FALSE ~ "Don't impute missing values (not required)",
    Missing_values==TRUE & Required==TRUE & Comorbidity==TRUE ~ "Don't impute missing values (comorbidity)",
    Missing_values==TRUE & Required==TRUE & Comorbidity==FALSE & Core==FALSE ~ "Don't impute missing values (no core variable)",
    Missing_values==TRUE & Required==TRUE & Comorbidity==FALSE & Core==TRUE & All_missing_values==TRUE ~ "Don't impute missing values (variable completely missing)",
    Missing_values==TRUE & Required==TRUE & Comorbidity==FALSE & Core==TRUE & All_missing_values==FALSE & cond_mcar==TRUE & Perc_miss_lt==TRUE ~ "Listwise deletion (MCAR)",
    Missing_values==TRUE & Required==TRUE & Comorbidity==FALSE & Core==TRUE & All_missing_values==FALSE & cond_mcar==TRUE & Perc_miss_lt==FALSE ~ "Imputation of missing values (MCAR)",
    Missing_values==TRUE & Required==TRUE & Comorbidity==FALSE & Core==TRUE & All_missing_values==FALSE & cond_mcar==FALSE & Perc_miss_lt5==TRUE ~ "Imputation of missing values (not MCAR)",
    Missing_values==TRUE & Required==TRUE & Comorbidity==FALSE & Core==TRUE & All_missing_values==FALSE & cond_mcar==FALSE & Perc_miss_lt5==FALSE & Perc_miss_lt15==TRUE ~ "Imputation of missing values (not MCAR, need for sensitivity analysis)",
    Missing_values==TRUE & Required==TRUE & Comorbidity==FALSE & Core==TRUE & All_missing_values==FALSE & cond_mcar==FALSE & Perc_miss_lt==FALSE & Perc_miss_lt15==FALSE ~ "Imputation of missing values (not MCAR, limit causal interpretation)"
  )) 

```

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

## 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 (MCAR)")),"Variable_name"])

## Description: variables that require listwise deletion
if(length(v_ld)>0) {
  tryCatch(
    {
      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=auxilary_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)
    },
    finally={
      dbDisconnect(con, shutdown=TRUE)
    })
}

```

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

## 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)","Imputation of missing values (not MCAR)","Imputation of missing values (not MCAR, need for sensitivity analysis)","Imputation of missing values (not MCAR, limit causal interpretation)")),"Variable_name"])

## Description: function for variables that require imputation
method_imputation <- function(df_inp, v) {
  tmp <- df_inp %>% select(-person_id) %>% mutate(residence_area_cd=as.factor(residence_area_cd))
  # Imputation age: if sex_cd!=2 or age_nm<12 or age_nm>50 -> impute pregnancy_bl to false
  tmp <- tmp %>% mutate(pregnancy_bl=ifelse(sex_cd!=2 | age_nm<12 | age_nm>50,FALSE,pregnancy_bl))
  if(length(v)>0) {
    ## If variables selected for imputation
    # -> Create objects
    pred <- make.predictorMatrix(tmp)
    meth <- make.method(tmp)
    meth[which(meth!="")] <- "pmm"
    # -> Run imputation
    imp <- mice(
      tmp,
      method = meth, 
      predictorMatrix = pred,
      m = 1, 
      maxit = 20
      )
    data_imp <- complete(imp) 
    data_imp <- cbind(person_id = df_inp$person_id, data_imp)
    # -> Select patient_ids of patients with imputed values
    out <- lapply(v, FUN =  function(variable) {
      out <- data.frame(person_id=df_inp[which(is.na(df_inp[,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_inp[0,]
    imp <- NULL
  }
  return(list(data_imp,imp))
}
# Use the function to obtain an imputed dataset
tmp <- method_imputation(df_not_viol %>% select(c(person_id, core_v)),v_imp)
data_imp <- tmp[[1]]

## Create database table with imputed data only
tryCatch(
    {
      con = dbConnect(duckdb::duckdb(), dbdir=auxilary_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"))
    },
    finally={
      dbDisconnect(con, shutdown=TRUE)
    })

```

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

## 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","Required","Comorbidity","Core","All_missing_values","MCAR","Perc_miss_lt","Perc_miss_lt5","Perc_miss_lt15")]

## 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(
              `Imputation method` = formatter("imputation_method",
                style = ~style(display = "block", 
                                padding = "0 4px", 
                                `border-radius` = "4px", 
                                `background-color` = case_when(`Imputation method`=="Imputation of missing values (MCAR)" ~"#F4A460",
                                                               `Imputation method`=="Imputation of missing values (not MCAR)" ~"#F4A460",
                                                               `Imputation method`=="Imputation of missing values (not MCAR, need for sensitivity analysis)" ~"#F4A460",
                                                               `Imputation method`=="Imputation of missing values (not MCAR, limit causal interpretation)" ~"#F4A460",
                                                               `Imputation method`=="No missing values" ~ "#E5FFCC",                                                                                                                                            `Imputation method`=="Don't impute missing values (variable completely missing)" ~ "#E0E0E0",
                                                               `Imputation method`=="Don't impute missing values (no core variable)" ~ "#E0E0E0",
                                                               `Imputation method`=="Don't impute missing values (comorbidity)" ~ "#E0E0E0",
                                                               `Imputation method`=="Don't impute missing values (not required)" ~ "#E0E0E0",
                                                               `Imputation method`=="Listwise deletion (MCAR)" ~ "#CCE5FF"
                                                               ))),
              `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("Imputation method" = imputation_method, "Number of imputed values" = Nr_imputed)


```

##### Distribution imputed data

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

imp <- tmp[[2]]

if(!is.null(imp)) {
  ## Create plots
  plot(imp)
  formula <- as.formula(paste("~", paste(v_imp, 
        collapse = "+", sep = ""), sep = ""))
  densityplot(imp,formula)
}

```
:::

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

tryCatch(
  {
    con = dbConnect(duckdb::duckdb(), dbdir=auxilary_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")

  },
  finally={
    dbDisconnect(con, shutdown=TRUE)
  })

formattable(count_ld,
            align="c")

```
