#####################################-
# Project: COVID Intubation Modeling - methods
# Analyst: Ryan Peterson
# Date Created: 2020-06-25
#####################################-

source("JDS_Submission/milasso.R")

library(tidyverse)
library(survival)
library(mice)
library(knitr)
library(xtable)
library(Hmisc)
library(survminer)
library(ncvreg)
library(pROC)
library(arsenal)
library(CIDAtools)

#### Data read/clean ####

df <- readxl::read_xlsx("DataRaw/Patient Data Set 4-9.xlsx", 1)
df$lddimer <- log(df$ddimer)
df$lferritin <- log(as.numeric(df$ferritin))

## properly label data 
df_mod <- df %>% 
  select(age, male:nosocomial, mp_gte_05, bmi:ldh, lddimer, crp, lferritin, intubated, neutrophils:nlr) %>% 
  mutate_at(vars(-age, -symptom_duration_ba, -bmi:-ldh, -crp, -lddimer, -lferritin, -neutrophils:-nlr), factor) 

##### Missing data imputation ##### 

## Imputing missing values 
## rf vals increase trees to 50

set.seed(199)
dfi <- mice(df_mod, m = 20, method = "rf", ntree = 50)
dfi

save(dfi, file = "DataProcessed/imputed.RData")

#### Predictive modeling ####

if(!exists("dfi"))
  load("DataProcessed/imputed.RData")

m <- dfi$m
f <- ~ scale(crp) + scale(ldh) + scale(lddimer) +
  scale(lferritin) + symptom_duration_ba + scale(age) + scale(bmi) +
  nonwhite + dm + male + scale(nlr)

# Data in format preffered by MI.lasso
dfi2 <- lapply(1:m, function(i) {
  df <- mice::complete(dfi, i)
  cbind(model.matrix(f, df)[,-1], as.numeric(as.character(df$intubated)))
})

### MI.lasso procedure ###

lamvec <- (2^(seq(-3,4,len=101)))^2/2

## select the optimal tuning parameter by BIC
bicvec = 0
for (i in 1:length(lamvec)) {
  fit = MI.lasso(dfi2, m, lambda=lamvec[i])
  bicvec[i] = fit$bic
}

## select tuning parameter with the smallest BIC
opt = order(bicvec)[1]

## fit the model with the optimal lambda
fit.opt = MI.lasso(dfi2, m, lambda=lamvec[opt])

### Perform MALCoM methods ##
set.seed(1)
b_mce <- b_mcm <- b_mcl <- 
  matrix(nrow = m, ncol = 12)
for(i in 1:m) {
  df_c <- complete(dfi,i)
  
  y <- as.numeric(as.character(df_c$intubated))
  X <- model.matrix(f, data = df_c)[, -1]
  
  fit <- cv.ncvreg(X,  y, family = 'binomial', penalty = "lasso")
  b_mcl[i,] <- coef(fit)
  
  fit <- cv.ncvreg(X,  y, family = 'binomial')
  b_mcm[i,] <- coef(fit)
  
  fit <- cv.ncvreg(X,  y, family = 'binomial', alpha = .25, 
                   penalty = "lasso")
  b_mce[i,] <- coef(fit)
}

malcom_b_mcm <- apply(b_mcm, 2, median)
malcom_b_mce <- apply(b_mce, 2, median)
malcom_b_mcl <- apply(b_mcl, 2, median)
b_mil <- apply(fit.opt$coefficients,1,mean)

names(malcom_b_mcm) <- names(malcom_b_mce) <- names(malcom_b_mcl) <- names(b_mil) <- 
  c("Intercept", colnames(X))

coefs <- cbind(malcom_b_mcl, b_mil, malcom_b_mcm, malcom_b_mce)


### Do linear approximate models for Supp table

### Perform MALCoM methods ##

set.seed(1)
b_app_mce <- b_app_mcm <- b_app_mcl <- 
  matrix(nrow = m, ncol = 12)
for(i in 1:m) {
  df_c <- complete(dfi,i)
  
  y <- as.numeric(as.character(df_c$intubated))
  X <- model.matrix(f, data = df_c)[, -1]
  
  fit <- cv.ncvreg(X,  y, penalty = "lasso")
  b_app_mcl[i,] <- coef(fit)
  
  fit <- cv.ncvreg(X,  y)
  b_app_mcm[i,] <- coef(fit)
  
  fit <- cv.ncvreg(X,  y,  alpha = .25, 
                   penalty = "lasso")
  b_app_mce[i,] <- coef(fit)
}

malcom_b_app_mcm <- apply(b_app_mcm, 2, median)
malcom_b_app_mce <- apply(b_app_mce, 2, median)
malcom_b_app_mcl <- apply(b_app_mcl, 2, median)

names(malcom_b_app_mcm) <- names(malcom_b_app_mce) <- names(malcom_b_app_mcl) <- names(b_mil) <- 
  c("Intercept", colnames(X))

coefs_app <- cbind(malcom_b_app_mcl, b_mil, malcom_b_app_mcm, malcom_b_app_mce)

##### Validation ##### 

## Get validation data 
df_mod2 <- df %>% 
  select(intubated, deceased, discharged, age, bmi,  male:nosocomial, 
         ldh, ddimer, crp, ferritin, neutrophils:nlr) %>% 
  mutate(Cohort = "Training")

## Get validation data 

valdf <- readxl::read_xlsx("DataRaw/Validation Cohort.xlsx", 2) 
valdf$days_to_intubation[is.na(valdf$days_to_intubation)] <- 5.5
valdf$days_to_intubation[valdf$days_to_intubation == 0] <- 0.1
valdf$days_to_intubation[is.na(valdf$days_to_intubation) & valdf$intubated == 1]  <- 0
valdf$deceased <- ifelse(is.na(valdf$deceased), 0, valdf$deceased)
valdf$lddimer <- log(valdf$ddimer)
valdf$lferritin <- log(valdf$ferritin)
valdf$Cohort <- "Validation"

df_full <- bind_rows(df_mod2, valdf) %>%  
  mutate_at(vars(-age, -bmi, -symptom_duration_ba, -ldh, -crp, -neutrophils:-nlr,
                 -ddimer, -ferritin, -intubated), factor) %>% 
  mutate_if(function(x) is.factor(x) & identical(levels(x), c("0",1)), 
            function(x) factor(x, 0:1, c("No", "Yes"))) %>%
  mutate(ethnicity = factor(ethnicity, 0:4, c("White", "Black", "Hispanic", "Asian", "Other"))) %>% 
  select(-id:-lferritin)

## Load dictionary
dict <- readxl::read_xlsx("DataRaw/Patient Data Set 4-9.xlsx", 2)
labs <- lapply(names(df_full), function(x) dict$Label[dict$Codename == x])
names(labs) <- names(df_full)
labs$Cohort <- "Cohort"
labs$intubated <- "Intubated"
labs$active_tobacco <- "Tobacco Use"
labs$nonwhite <- "Nonwhite"
labs$nlr <- "NLR"

tab1 <- tableby(Cohort ~ ., data = df_full, 
                numeric.simplify = T, numeric.test = "kwt", 
                cat.test = "fe", digits = 2,
                simulate.p.value = TRUE,
                cat.simplify = T,  numeric.stats = c("meansd"), 
                cat.stats = c("countpct"), total = FALSE, test = F)
tab1_pretty <- as.data.frame(summary(tab1, labelTranslations = labs, text = "latex"))

# tab1_pretty %>%
#   kable(escape = F, align = "lrr") %>%
#   kable_styling(c("striped", "condensed"), full_width = F, fixed_thead = T) %>%
#   column_spec(1, width = "3cm")

## Validating predictive model 

# load("DataProcessed/mcp-fits.RData")

new_df <- valdf %>% 
  select(intubated, crp, ldh, lddimer, lferritin, symptom_duration_ba, age, bmi, nonwhite, dm, male, nlr) 

## Should be no missings in CRP, LDH, DM
new_df %>% 
  apply(2, function(x) sum(is.na(x)))

# Use mean imputation
new_df <- new_df %>% 
  mice::mice(m=1, method = "mean") %>% 
  complete()

y <- as.numeric(new_df$intubated)

new_X <- model.matrix(f, new_df)[,-1]

head(new_X)

lodds <- c(cbind(1, new_X) %*% malcom_b_mcl)
probs_mcl <- exp(lodds) / (1 + exp(lodds))

lodds <- c(cbind(1, new_X) %*% malcom_b_mcm)
probs_mcm <- exp(lodds) / (1 + exp(lodds))

lodds <- c(cbind(1, new_X) %*% malcom_b_mce)
probs_mce <- exp(lodds) / (1 + exp(lodds))

lodds <- c(cbind(1, new_X) %*% b_mcl[1,])
probs_scl <- exp(lodds) / (1 + exp(lodds))

probs_mil <- c(cbind(1, new_X) %*% b_mil)

g <- pROC::roc(y ~ probs_mcl, percent = T, ci = T, quiet = T)
g2 <- pROC::roc(y ~ probs_mil, percent = T, ci = T, quiet = T)
g3 <- pROC::roc(y ~ probs_mcm, percent = T, ci = T, quiet = T)
g4 <- pROC::roc(y ~ probs_mce, percent = T, ci = T, quiet = T)
g5 <- pROC::roc(y ~ probs_scl, percent = T, ci = T, quiet = T)

par(pty="s", mfrow = c(1,1))

plot(g, print.auc  = T)    
plot(g2, add = T, col = 2, print.auc = T, print.auc.adj = c(0,2.5))
plot(g3, add = T, col = 3, print.auc = T, print.auc.adj = c(0,4))
plot(g4, add = T, col = 4, print.auc = T, print.auc.adj = c(0,5.5))
# plot(g5, add = T, col = 4, print.auc = T, print.auc.adj = c(0,7.0))

legend("bottomright", c("MCL", "MIL", "MCM", "MCE"), col = 1:4, lwd = 2)

#### Likelihood 

my_dev <- function(y, pi_hat) 
  -2*sum(y * log(pi_hat) + (1-y) * log(1-pi_hat))

D0 <- my_dev(y, mean(valdf$intubated))

D1 <- list(
  mcl = my_dev(y, probs_mcl),
  mil = my_dev(y, probs_mil),
  mcm = my_dev(y, probs_mcm),
  mce = my_dev(y, probs_mce)
)
  
R2 <- 100*sapply(D1, function(D) (D0- D)/D0)

Accuracy <- 100*c(
  mean((probs_mcl >= .5) == (y == 1)), 
  mean((probs_mil >= .5) == (y == 1)), 
  mean((probs_mcm >= .5) == (y == 1)), 
  mean((probs_mce >= .5) == (y == 1))
)

AUC <- c(
  g$auc, 
  g2$auc, 
  g3$auc, 
  g4$auc
)

(accuracy_table <- cbind(AUC, R2, Accuracy))

## Summary stats 


df <- readxl::read_xlsx("DataRaw/Patient Data Set 4-9.xlsx", 1)
df$days_to_intubation[!is.na(df$days_to_intubation_indicated)] <- df$days_to_intubation_indicated[!is.na(df$days_to_intubation_indicated)]
df$days_to_intubation[is.na(df$days_to_intubation)] <- 5.5
df$days_to_intubation[df$days_to_intubation == 0] <- 0.1

## From meeting, reset those with missing max pressors to zero
df$mp_gte_05 <- ifelse(is.na(df$mp_gte_05), 0, df$mp_gte_05)

df$days_to_intubation[is.na(df$days_to_intubation) & df$intubated == 1]  <- 0
df$deceased <- ifelse(is.na(df$deceased), 0, df$deceased)

# From email, for patients who died before intubation, enter in date when intubation first indicated
df$days_to_intubation[!is.na(df$days_to_intubation_indicated)] <- df$days_to_intubation_indicated[!is.na(df$days_to_intubation_indicated)]

## For survival analysis:
df$days_to_intubation[is.na(df$days_to_intubation)] <- 5.5
df$days_to_intubation[df$days_to_intubation == 0] <- 0.1
df$lddimer <- log(df$ddimer)
df$lferritin <- log(df$ferritin)

df$intubated <- factor(df$intubated, 0:1, c("Not intubated", "Intubated"))

## properly label data 
df_mod <- df %>% 
  select(intubated, deceased, discharged, age, bmi,  male:nosocomial, ldh, ddimer, crp, ferritin, neutrophils:nlr) %>% 
  mutate(Cohort = "Training")


## Get validation data 

valdf <- readxl::read_xlsx("DataRaw/Validation Cohort.xlsx", 2) 
valdf$days_to_intubation[is.na(valdf$days_to_intubation)] <- 5.5
valdf$days_to_intubation[valdf$days_to_intubation == 0] <- 0.1
valdf$days_to_intubation[is.na(valdf$days_to_intubation) & valdf$intubated == 1]  <- 0
valdf$deceased <- ifelse(is.na(valdf$deceased), 0, valdf$deceased)
valdf$lddimer <- log(valdf$ddimer)
valdf$lferritin <- log(valdf$ferritin)
valdf$intubated <- factor(valdf$intubated, 0:1, c("Not intubated", "Intubated"))
valdf$Cohort <- "Validation"

df_full <- bind_rows(df_mod, valdf) %>%  
  mutate_at(vars(-age, -bmi, -symptom_duration_ba, -ldh, -crp, -neutrophils:-nlr, -ddimer, -ferritin, -intubated), factor) %>% 
  mutate_if(function(x) is.factor(x) & identical(levels(x), c("0",1)), function(x) factor(x, 0:1, c("No", "Yes"))) %>%
  mutate(ethnicity = factor(ethnicity, 0:4, c("White", "Black", "Hispanic", "Asian", "Other"))) %>% 
  select(-id:-lferritin, -deceased, -discharged, -cr_gt2, -starts_with("cx_"), -ever_proned, -ever_pressors, -nosocomial)

tab1 <- tableby(Cohort ~ ., data = df_full, 
                numeric.simplify = T, numeric.test = "kwt", 
                cat.test = "fe", digits = 2,
                simulate.p.value = TRUE,
                cat.simplify = T,  numeric.stats = c("meansd"), 
                cat.stats = c("countpct"), total = FALSE, test = F)
tab1_pretty <- as.data.frame(summary(tab1, labelTranslations = labs, text = "latex"))

tab1_pretty %>% 
  kable(escape = F, align = "lrr") %>% 
  kable_styling(c("striped", "condensed"), full_width = F, fixed_thead = T) %>% 
  column_spec(1, width = "3cm") 

## Table 2: Missing data by intubation status

t2_data <- df_mod %>% 
  select(crp, ldh, ddimer, ferritin, symptom_duration_ba, age, bmi, nonwhite, dm, male, nlr) %>% 
  mutate_all(function(x) factor(1*is.na(x), 0:1, c("Nonmissing", "Missing"))) %>% 
  mutate("Any of the above" = apply(., 1, function(x) any(x == "Missing"))) %>%
  bind_cols(intubated = df_mod$intubated) 
  
tab2 <- tableby(intubated ~ ., data = t2_data, 
                numeric.simplify = T, numeric.test = "kwt", 
                cat.test = "fe", digits = 2,
                simulate.p.value = TRUE,
                cat.simplify = T,  numeric.stats = c("Nmiss", "meansd", "medianq1q3"), 
                cat.stats = c("Nmiss", "countpct"), total = T)
tab2_pretty <- as.data.frame(summary(tab2, labelTranslations = labs, text = "latex"))

# tab2_pretty %>% 
#   kable(escape = F, align = "lrr") %>% 
#   kable_styling(c("striped", "condensed"), full_width = F, fixed_thead = T) %>% 
#   column_spec(1, width = "3cm") 

### Table 4 with comparable values



save(f, tab1_pretty, tab2_pretty, coefs, b_mcm, b_mce, b_mcl, coefs_app,
     accuracy_table, g, g2, g3, g4, file= "JDS_Submission/application.RData")
