misc.Rmd

---
title: "Investigating the epidemiology of MISC"
output:
  word_document:
    toc: yes
  html_document:
    code_folding: hide
    toc: yes
    toc_float: yes
editor_options:
  markdown:
    wrap: 72
---

# Setting up the R environment

### Installing packages and loading the library

```{r libraries, message=FALSE, warning=FALSE, results='hide'}
# Install packages
paket <- function(pak){
  new_pak <- pak[!(pak %in% rownames(installed.packages()))]
  if (length(new_pak)) 
    install.packages(new_pak, dependencies = TRUE,repos="https://cloud.r-project.org/")
  sapply(pak, library, character.only = TRUE)
}

listOfPackages <- c("tidyverse", "RColorBrewer", "knitr", "kableExtra","tableone", "gridExtra", "dplyr", "lubridate", "MatchIt", "table1")
paket(listOfPackages)
```

### R session information

```{r}
sessionInfo()
```

### Plots aesthetics information

```{r}
theme_plots <- theme_bw() +
  theme(strip.text = element_text(size = 5),
        axis.text.x = element_text(size = 8), 
        axis.text.y = element_text(size = 6), 
        axis.title.x = element_text(size = 8),
        axis.title.y = element_text(size = 8),
        title = element_text(size = 10),
        plot.subtitle = element_text(size = 9, face = "italic")) 
theme_set(theme_plots)

# Colorblind palette
cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
```

# Data retrieval

## Inclusion criteria
Patients inclusion criteria:

    -   Children (<21 years) 
    -   Diagnosed with MIS-C based on CDC/WHO/RCPCH criteria  


## Variables that need to be checked/modified by each site

Change the values of the following variables according to the specificities of your site:

1.  "folder_4ce_files": folder path where your phase 2.2 data files are
    located

2.   "obfuscation": determine the obfuscation threshold (FALSE if no
    obfuscation, or the numeric value of the obfuscation threshold if
    any)

3.   "raceAvailable": set as TRUE or FALSE depending on whether the
    variable is being collected at your site

4.   "dateFormat": specify the format of the date at your site (e.g., for
    "03-AUG-20", the format would be "%d-%b-%y", [see
    documentation](https://www.stat.berkeley.edu/~s133/dates.html))
    
5. "data_update_date": date at which the data has been updated in the local 
data warehouse. Used to estimate patient age at time of visit, since patients age
in the 4CE demographic file is expected the age at data update. 

6. "country": 

```{r message=FALSE, warning=FALSE}
folder_4ce_files <- "/4ceData/Input/4ce/"
obfuscation =  FALSE
raceAvailable = TRUE
dateFormat <- "%d-%b-%y"
data_update_date <- "2022-04-02"
country <- "US"
```

## Data loading

We will use as input the 2.2 data files. Specifically: 

- LocalPatientSummary 

- LocalPatientObservation 

- LocalPatientClinicalcourse

For sites recording the race, we will also use an additional file: 

- LocalPatientRace

### Reading 4CE phase 2.2 files

```{r message=FALSE, warning=FALSE}
### Read the CSV input files
source("R/readInputFiles.R")

getwd()
files <- readInputFiles( path      = folder_4ce_files, 
                         separator = ",",
                         skip      = 0, 
                         verbose   = FALSE )
  
## Create the output folder if it doesn't exist
if (! "output" %in% list.dirs()) dir.create("output")

### Extract the patient summary and observation information. 
demo_raw <- files[["patientSummary"]] 
obs_raw <- files[["patientObservations"]] 
clinical_raw <- files[["patientClinicalCourse"]] 

### Read the file containing race information for those sites recording this variable
if( raceAvailable == TRUE ){
  race_raw <- read.delim(file.path(folder_4ce_files, "/LocalPatientRace.csv"), sep = ",", skip = 0)
}

### Read the file containing the variants dates per country
variantsDates <- read.delim("./public-data/variantsDates.txt", sep = "\t") %>%
  filter( Country == country )
```

## Extract the MIS-C patients and assign the variant based on the date they were diagnosed with MISC
```{r}
misc_patients <- obs_raw %>%
  filter( concept_type == "COVID-MISC") 

misc_all <- left_join( misc_patients, clinical_raw[, c("patient_num", "days_since_admission", "calendar_date", "in_hospital", "severe", "in_icu", "dead")], by=c("patient_num", "days_since_admission"))

misc_all <- left_join( misc_all, demo_raw[, c("patient_num", "age", "sex", "admission_date")], by= "patient_num") %>%
  mutate( date = as.Date( admission_date, format = dateFormat ), 
          variant_misc = ifelse( date >= variantsDates$Omicron, "Omicron", ifelse( date <= variantsDates$Alpha, "Alpha", "Delta"))) 

misc_complete <- misc_all %>%
  filter( !is.na( calendar_date )) %>%
   dplyr::mutate( weeks = as.Date(cut( date, breaks = "week")),
                  month = as.Date(cut( date, breaks = "month")),
                  year = format( date, "%Y"))
```

## Estimate number of MIS-C cases per month
```{r}
misc_cases <- misc_complete %>%
  dplyr::group_by( month ) %>%
  dplyr::summarise( distinct_patients = length(unique(patient_num))) 

### plot barplot
ggplot(data=misc_cases, aes(x=month, y=distinct_patients)) +
  geom_bar(stat="identity") 

## number of patients by misc variant
misc_complete %>%
  dplyr::group_by( variant_misc ) %>%
  dplyr::summarise( distinct_patients = length(unique(patient_num))) 
```

## Estimate for the MIS-C patients death and ICU admissions per wave 
```{r}
misc_icu <- misc_complete %>%
  filter( in_icu == 1) %>%
  dplyr::group_by( variant_misc ) %>%
  dplyr::summarise( icu_patients = length(unique(patient_num))) %>%
  dplyr::mutate( var = "ICU")

ggplot(misc_icu, aes(fill=variant_misc, y=icu_patients, x=var)) + 
    geom_bar(position="dodge", stat="identity") +
  scale_fill_manual( values = cbPalette)
```
## Estimate for the MIS-C patients age and sex ratio
```{r}
misc_complete %>%
  dplyr::group_by( variant_misc ) %>%
  ggplot(aes(variant_misc, age))+
  geom_boxplot(aes(fill=variant_misc), position = "dodge") +
  scale_fill_manual( values = cbPalette)


misc_sex_ratio <- misc_complete %>%
  dplyr::group_by( variant_misc, sex ) %>%
  dplyr::summarise( n =  length(unique(patient_num))) %>%
  dplyr::group_by( variant_misc ) %>%
  dplyr::mutate( ratio = n/sum(n))

ggplot(misc_sex_ratio, aes(fill=variant_misc, y=n, x=sex)) + 
    geom_bar(position="dodge", stat="identity") +
  scale_fill_manual( values = cbPalette)
```

## Estimate the length of the MISC hospitalization 

```{r}
### this code is copy-paste from the peds psy study
clinical_raw_hosp <- filter(clinical_raw, in_hospital == 1)

count_sequences_hospitalisation <- function(df, ...) {
  seq_hospitalisation_df <- data.frame(total_span = seq(min(df$days_since_admission),
                                                        max(df$days_since_admission))
  ) %>%
    left_join(df, by = c("total_span" = "days_since_admission")) %>%
    replace_na(list(in_hospital = 0))
  count_sequences <- rle(seq_hospitalisation_df$in_hospital)
  count_sequences_1 <- lapply(count_sequences, function(x) x[count_sequences$values == 1])
  n_sequences <- seq_along(count_sequences_1$lengths)
  sequences <- rep.int(n_sequences, count_sequences_1$lengths)
  sequences_len <- rep.int(count_sequences_1$lengths, count_sequences_1$lengths)
  stopifnot(length(df$days_since_admission) == length(sequences))
  data.frame(days_since_admission = df$days_since_admission,
             n_hospitalisation = sequences,
             len_hospitalisation = sequences_len)
}
stopifnot(all(clinical_raw_hosp$in_hospital == 1))
hospitalisations_seq_df <- clinical_raw_hosp %>%
  distinct(patient_num, cohort, days_since_admission, in_hospital) %>%
  group_by(patient_num, cohort) %>%
  group_modify(count_sequences_hospitalisation)

misc_hosp_length <- left_join(misc_complete,
                          hospitalisations_seq_df,
                          by = c("patient_num", "cohort", "days_since_admission"))

misc_hosp_length$len_hospitalisation <- ifelse( is.na(misc_hosp_length$len_hospitalisation), 0, misc_hosp_length$len_hospitalisation)
misc_hosp_length %>%
  dplyr::group_by( variant_misc ) %>%
  ggplot(aes(variant_misc, as.numeric(len_hospitalisation)))+
  geom_boxplot(aes(fill=variant_misc), position = "dodge", outlier.shape = NA) +
  scale_y_continuous(limits = quantile(misc_hosp_length$len_hospitalisation,na.rm = T,  c(0.1, 0.9))) +
  scale_fill_manual( values = cbPalette)

```

### Summary of procedures and meds present in the MISC patients during MISC hospitalization
```{r}
misc_hosp_length <- misc_hosp_length %>%
  mutate( comb = paste0( cohort, patient_num, n_hospitalisation ) )

misc_all_days_hospitalized <- hospitalisations_seq_df %>%
    mutate( key = paste0( cohort, patient_num, days_since_admission ), 
            comb = paste0( cohort, patient_num, n_hospitalisation ) ) %>%
  filter( comb %in% misc_hosp_length$comb)


## select all days of the hospitalization

obs_data_filter <- obs_raw %>%
  filter( patient_num %in% misc_complete$patient_num ) %>%
  mutate( key = paste0( cohort, patient_num, days_since_admission )) %>%
  filter( key %in% misc_all_days_hospitalized$key )

meds_sum <- obs_data_filter %>% 
  filter( concept_type == "MED-CLASS") %>% 
  dplyr::group_by( concept_code ) %>% 
  dplyr::summarise( patients =  length(unique(patient_num))) %>%
  dplyr::arrange( desc(patients) )

print( meds_sum) 

prod_sum <- obs_data_filter %>% 
  filter( concept_type == "PROC-GROUP") %>% 
  dplyr::group_by( concept_code ) %>%
  dplyr::summarise( patients =  length(unique(patient_num))) %>%
  dplyr::arrange( desc(patients) )

print( prod_sum)
```



## Prepare a table to perform logistic regression

```{r}
mainTable <- misc_hosp_length %>%
  select( patient_num, date, variant_misc, len_hospitalisation, age, sex, in_icu, dead )
        
if( raceAvailable == TRUE ){
  mainTable <- left_join( mainTable, race_raw %>% select( patient_num, race_4ce), by="patient_num")
}          

### add other variables present in the clinical characteristic table
clinicalCharact <- read.delim('public-data/clinicalCharacteristics.txt')

clin_var <- obs_data_filter %>%
  filter( concept_code %in% clinicalCharact$concept_code ) %>%
  left_join( clinicalCharact, by = "concept_code") %>%
  mutate( value = 1 ) %>%
  select( patient_num, concept_code, variableName, value) %>%
  unique()

clin_var2merge <- clin_var %>%
  spread(key = variableName, value =  value, fill = 0) %>%
  select( -concept_code)

### add the new variables and  
mainTable <- mainTable %>%
  dplyr::left_join( clin_var2merge, by = "patient_num") %>%
  replace(is.na(.), 0) 

## any non-present variable as empty column all with 0
varsNoPatients <- clinicalCharact %>%
  dplyr::select( variableName ) %>%
  unique( ) %>%
  dplyr::filter(! variableName %in% colnames( mainTable )) 

for( i in 1:length( varsNoPatients$variableName )){
  mainTable[ varsNoPatients$variableName[i]] <- 0
}

## Add a column called cardiovascular outcome
# logic to fill this new colum: it will be a 1 if the patient has 
# Cardiovascular outcome: SICARDIAC OR/AND Shock related ICD-10 code AND/OR ECMO  AND/OR CARDIAC ARREST AND/OR CPR
mainTable <- mainTable %>%
  dplyr::mutate( cardiovascular_outcome = ifelse( SICARDIAC == 1 | Shock == 1| ECMO == 1 | `Cardiac arrest` == 1 | CPR == 1, 1, 0))

## check if all patients with cardiovascular outcome where in ICU
# toCheck <- mainTable %>% 
#  dplyr::filter( cardiovascular_outcome == 1) %>%
#  dplyr::mutate( check = ifelse( in_icu == 1, "OK", "CHECK" )) %>%
#  dplyr::filter( check == "CHECK") %>%
#  dplyr::select( patient_num, in_icu, cardiovascular_outcome, SICARDIAC, Shock,ECMO, CPR, `Cardiac arrest`  )
```

## Table 1
```{r}
mainTable$sex    <- factor(mainTable$sex, levels=c("male", "female"), labels=c("Male", "Female"))
mainTable$variant_misc    <- factor(mainTable$variant_misc)
mainTable$race_4ce     <- factor(mainTable$race_4ce, 
                                 levels=c("white", "black", "asian", "american_indian", "other", "n_information"),
                                 labels=c("White", "Black", "Asian", "American indian", "Other", "Not Available"))

mainTable$in_icu  <- as.logical(mainTable$in_icu == 1)
mainTable$dead <- as.logical(mainTable$dead == 1)
mainTable$Shock <- as.logical(mainTable$Shock == 1)
mainTable$ECMO <- as.logical(mainTable$ECMO == 1)
mainTable$SICARDIAC <- as.logical(mainTable$SICARDIAC == 1)
mainTable$CPR <- as.logical(mainTable$CPR == 1)
mainTable$cardiovascular_outcome <- as.logical(mainTable$cardiovascular_outcome == 1)
mainTable$`Heart failure` <- as.logical(mainTable$`Heart failure` == 1)
mainTable$`Respiratory failure` <- as.logical(mainTable$`Respiratory failure` == 1)
mainTable$`Cardiac arrest` <- as.logical(mainTable$`Cardiac arrest` == 1)

label(mainTable$race_4ce) <- "Race"
label(mainTable$age)      <- "Age"
label(mainTable$len_hospitalisation) <- "Length of hospitalization (days)"
label(mainTable$sex)  <- "Sex"
label(mainTable$in_icu) <- "ICU"
label(mainTable$cardiovascular_outcome)     <- "Cardiovascular Outcome"
label(mainTable$variant_misc)     <- "Variant MIS-C"

```

## Function to estimate the p-values extracted from the table1 vignettes
```{r}
# using as reference https://cran.r-project.org/web/packages/table1/vignettes/table1-examples.html 
ownpvalue <- function(x, ...) {
    # Construct vectors of data y, and groups (strata) g
    y <- unlist(x)
    g <- factor(rep(1:length(x), times=sapply(x, length)))
    if (is.numeric(y)) {
        # For numeric variables, perform a standard 2-sample t-test
        p <- t.test(y ~ g)$p.value
    } else {
        # For categorical variables, perform a chi-squared test of independence
        p <- chisq.test(table(y, g))$p.value
    }
    # Format the p-value, using an HTML entity for the less-than sign.
    # The initial empty string places the output on the line below the variable label.
    c("", sub("<", "&lt;", format.pval(p, digits=3, eps=0.001)))
}
```

## Create the table
```{r}
table1::table1(~ age + race_4ce + len_hospitalisation + in_icu + dead + cardiovascular_outcome + Shock + SICARDIAC + ECMO + `Heart failure` + `Cardiac arrest` + `Respiratory failure` | variant_misc,
    data=mainTable, extra.col=list(`P-value`=ownpvalue))
```

### Check if we have mechanical ventilation or intubation 

```{r}
# procedures_all <- obs_raw %>% 
#   filter( concept_type == "PROC-GROUP") %>% 
#   dplyr::group_by( concept_code ) %>%
#   dplyr::summarise( patients =  length(unique(patient_num))) %>%
#   dplyr::arrange( desc(patients) )
# 
# obs_data_filter %>% filter( concept_code == "Z99.11")
```

### Laboratory values during the MIS-C hospitalization

```{r}
### read the table with the labs of interest and the expected ranges
laboratoryCharact <- read.delim('public-data/laboratoryCharacteristics.txt')

lab_var <- obs_data_filter %>%
  filter( concept_code %in% laboratoryCharact$concept_code ) %>%
  left_join( laboratoryCharact, by = "concept_code") %>%
  select( patient_num, concept_code, variableName, value, days_since_admission )

getMaxValues <- lab_var %>%
  dplyr::group_by( patient_num, concept_code ) %>%
  dplyr::mutate( max = max( value ), 
                 min = min( value ))
```