GLM_MCL.Rmd

---
title: "GLM Model"
output: html_notebook
author: HA Hill
---


```{r}
library(tidymodels)

```


#Make a Multivariate Logistic Model with Top 20 variable importance
factors from 10 fold CV (full model).

##Split into train and test set - use 10-cross-validation for parameter
fit

```{r}
mcl_glm <- mcl_glm %>%  mutate_if(is.character, as.factor)
mcl_glm$status <- relevel(mcl_glm$status, ref = "mild")

set.seed(102)
mcl_glm_split <- initial_split(mcl_glm, strata = status)
mcl_glm_train <- training(mcl_glm_split)
mcl_glm_test <- testing(mcl_glm_split)
mcl_glm_folds <- vfold_cv(mcl_glm_train, strata=status)
```

#Will need to impute all predictor variables - use k nearest neighbor
imputation (KNN).
#Normalize all numeric predictors

```{r}
mcl_glm_recipe <- recipe(status ~., data=mcl_glm_train) %>% 
  update_role(pt_id, new_role = "id variable" ) %>% 
  step_other(c(morph, ecog), threshold = 0.01 ) %>% 
  step_impute_knn(all_predictors()) %>% 
  step_normalize(all_numeric_predictors()) %>% 
  step_dummy(all_nominal_predictors()) 
  

 

mcl_glm_recipe_prep <- prep(mcl_glm_recipe)
juice(mcl_glm_recipe_prep)

```
#Set computational engine.

```{r}
mcl_glm_spec <- logistic_reg() %>% 
  set_engine("glm")

mcl_glm_spec
```

#Put together workflow

```{r}
mcl_glm_wf <- workflow() %>% 
  add_recipe(mcl_glm_recipe)

mcl_glm_wf
```

#Fit best model on validation and training sets.

```{r}
mcl_glm_metrics <- metric_set(roc_auc, accuracy)

doParallel::registerDoParallel()

mcl_glm_rs <- mcl_glm_wf %>% 
  add_model(mcl_glm_spec) %>% 
  fit_resamples(
    resamples=mcl_glm_folds,
    metrics = mcl_glm_metrics,
    control=control_resamples(save_pred = TRUE)
  )

mcl_glm_rs
```

#Fit best model on test set. 
##Collect metrics, predictions and visualize confustion
matrix.

```{r}
collect_metrics(mcl_glm_rs)
```

```{r}
mcl_glm_final <- mcl_glm_wf %>% 
  add_model(mcl_glm_spec) %>% 
  last_fit(mcl_glm_split)

mcl_glm_final
collect_metrics(mcl_glm_final)
```

```{r}
collect_predictions(mcl_glm_final) %>% 
  conf_mat(status, .pred_class) %>% 
  autoplot()
```

```{r}
collect_metrics(mcl_glm_final)
collect_predictions(mcl_glm_final) %>% 
  conf_mat(status, .pred_class)
```
#Get coefficients from logistic model

```{r}
mcl_glm_final$.workflow[[1]] %>% 
  tidy(exponentiate=TRUE)
```


```{r}
library(kableExtra)
library(gtsummary)
library(gt)

 glm_20 <- mcl_glm_rs %>%
  pull(.workflow) %>%
  pluck(1) %>%
  tidy() %>% 
  arrange(p.value) %>%
   gt() %>% 
   tab_header(
     title = "Multivariate Logistic Model") %>% 
fmt_number(columns = c(estimate, std.error, p.value, statistic),  decimals=2) %>%
    tab_style(
    style = list(
      cell_fill(color = "lightcyan"),
      cell_text(weight = "bold")
      ),
    locations = cells_body(
      columns = p.value,
      rows = p.value<= 0.06
    )
  ) 
       
```
#Make table with coefficients and results
```{r}
GLM_results %>%
  gt() %>%
  tab_header(
     title = "Multivariate Logistic Model") %>% 
  fmt_number(columns = c("Odds Ratio"),  decimals=2) %>% 
    tab_style(
    style = list(
      cell_fill(color = "lightcyan"),
      cell_text(weight = "bold")
      ),
    locations = cells_body(
      columns = "p-value",
      rows = "p-value"<= 0.06
    )
  ) 
```


```{r}

 glm_20
 
```

```{r}
cols_label(ki_67 = "ki67 %",
              ldh = "LDH",
              age_dx_years = "age at diagnosis",
              morph_leukemic = "leukemic morphology",
              total_genes_mutated = "total somatic mutations",
              smoke_yes="current smoker",
              bm = "bone marrow involvement",
              b2m = "beta 2 microglobulin",
              wbc = "wbc count",
              morph_blastoid="blastoid morpology",
              b_symptoms_yes="b symptoms",
              hgb="hemoglobin",
              notch1_yes = "NOTCH1 mutation",
              platelets="platelet count",
              smoke_no="never smoker",
              gi_involve_yes="GI involvement",
              birc3_yes="BIRC3 mutation",
              ros1_yes="ROS1 mutation",
              ecog_X2= "ECOG Performance Status 2",
              morph_blastoid_and_pleomorphic = "blastoid and pleomorphic      morphology",
              ecog_X1="ECOG Performance Status 1",
              morph_classic= "classic morphology",
              morph_pleomorphic="pleomorphic morphology",
              bmi = "BMI",
              card11_yes = "CARD11 mutation",
              ecog_other = "ECOG Performance Status 3 or 4",
              tp53_yes = "TP53 mutation"
           
              )

cols_label("ki_67" = "ki67 %",
              "ldh" = "LDH",
              "age_dx_years" = "age at diagnosis",
              "morph_leukemic" = "leukemic morphology",
              "total_genes_mutated" = "total somatic mutations",
              "smoke_yes"="current smoker",
              "bm" = "bone marrow involvement",
              "b2m" = "beta 2 microglobulin",
              "wbc" = "wbc count",
              "morph_blastoid"="blastoid morpology",
              "b_symptoms_yes"="b symptoms",
              "hgb"="hemoglobin",
              "notch1_yes" = "NOTCH1 mutation",
              "platelets"="platelet count",
              "smoke_no"="never smoker",
              "gi_involve_yes"="GI involvement",
              "birc3_yes"="BIRC3 mutation",
              "ros1_yes"="ROS1 mutation",
              "ecog_X2"= "ECOG Performance Status 2",
              "morph_blastoid_and_pleomorphic" = "blastoid and pleomorphic      morphology",
              "ecog_X1"="ECOG Performance Status 1",
              "morph_classic"= "classic morphology",
              "morph_pleomorphic"="pleomorphic morphology",
              "bmi" = "BMI",
              "card11_yes" = "CARD11 mutation",
              "ecog_other" = "ECOG Performance Status 3 or 4",
              "tp53_yes" = "TP53 mutation")
```


```{r}
tbl_regression(mcl_glm_rs)
   
```

```{r}
gtsave(glm_20, "GLM_20.docx")
```

#Plot model coefficients

```{r}
library(forcats)

mcl_glm_final %>%
  pull(.workflow) %>%
  pluck(1) %>%
  tidy() %>%
  filter(term != "(Intercept)") %>%
  ggplot(aes(estimate, fct_reorder(term, estimate))) +
  geom_vline(xintercept = 0, color = "black", lty = 2, size = 1.2) +
  geom_errorbar(aes(
    xmin = estimate - std.error,
    xmax = estimate + std.error
  ),
  width = .2, color = "black", alpha = 0.7
  ) +
  geom_point(size = 2, color = "blue") +
  labs(y = NULL, x = "Coefficient from logistic regression")
```

```{r}
glm_predict <- mcl_glm_final %>%
  collect_predictions() %>%
  group_by(id)
```

```{r}
write_xlsx(glm_predict, "glm_pred.xlsx")
```

#GLM model with 10 highest variables from full model.
##Split into train and test set - use 10-cross-validation for parameter
fit

```{r}
mcl_glm_10 <- mcl_glm_10 %>%  mutate_if(is.character, as.factor)

set.seed(102)
mcl_glm_split <- initial_split(mcl_glm_10, strata = status)
mcl_glm_train <- training(mcl_glm_split)
mcl_glm_test <- testing(mcl_glm_split)
mcl_glm_folds <- vfold_cv(mcl_glm_train, strata=status)
```

#Will need to impute all predictor variables - use k nearest neighbor
imputation (KNN).
#Normalize all numeric predictors

```{r}
mcl_glm_recipe <- recipe(status ~., data=mcl_glm_train) %>% 
  update_role(pt_id, new_role = "id variable" ) %>% 
  step_other(morph, ecog) %>% 
  step_impute_knn(all_predictors()) %>% 
  step_normalize(all_numeric_predictors()) %>% 
  step_dummy(all_nominal_predictors()) 
  

 

mcl_glm_recipe_prep <- prep(mcl_glm_recipe)
juice(mcl_glm_recipe_prep)

```
#Set computational engine.

```{r}
mcl_glm_spec <- logistic_reg() %>% 
  set_engine("glm")

mcl_glm_spec
```

#Put together workflow

```{r}
mcl_glm_wf <- workflow() %>% 
  add_recipe(mcl_glm_recipe)

mcl_glm_wf
```
#Fit best model on validation and training sets.

```{r}
mcl_glm_metrics <- metric_set(roc_auc, accuracy)

doParallel::registerDoParallel()

mcl_glm_rs <- mcl_glm_wf %>% 
  add_model(mcl_glm_spec) %>% 
  fit_resamples(
    resamples=mcl_glm_folds,
    metrics = mcl_glm_metrics,
    control=control_resamples(save_pred = TRUE)
  )

mcl_glm_rs
```
#Fit best model on test set. ##Collect metrics and visualize confustion
matrix.

```{r}
collect_metrics(mcl_glm_rs)
```

```{r}
mcl_glm_final <- mcl_glm_wf %>% 
  add_model(mcl_glm_spec) %>% 
  last_fit(mcl_glm_split)

mcl_glm_final
collect_metrics(mcl_glm_final)
```

```{r}
library(kableExtra)
library(gt)

 glm_20 <- mcl_glm_final %>%
  pull(.workflow) %>%
  pluck(1) %>%
  tidy() %>% 
  arrange(estimate) %>%
   gt() %>% 
   tab_header(
     title = "Coefficients from Multivariate Logistic Model"
   ) %>% 
   fmt_number(columns = c(estimate, std.error, p.value, statistic),  decimals=2)
 
gtsave(glm_20, "GLM_10.docx")
```

```{r}
collect_predictions(mcl_glm_final) %>% 
  conf_mat(status, .pred_class) %>% 
  autoplot()
```
```{r}
library(forcats)

mcl_glm_final %>%
  pull(.workflow) %>%
  pluck(1) %>%
  tidy() %>%
  filter(term != "(Intercept)") %>%
  ggplot(aes(estimate, fct_reorder(term, estimate))) +
  geom_vline(xintercept = 0, color = "black", lty = 2, size = 1.2) +
  geom_errorbar(aes(
    xmin = estimate - std.error,
    xmax = estimate + std.error
  ),
  width = .2, color = "black", alpha = 0.7
  ) +
  geom_point(size = 2, color = "#85144B") +
  labs(y = NULL, x = "Coefficient from logistic regression")`
```

```{r}
glm_predict <- mcl_glm_final %>%
  collect_predictions() %>%
  group_by(id)
```

```{r}
write_xlsx(glm_predict, "glm_pred_10.xlsx")
```

#Univariate logistic models with MIPIs
#Univariate models incorporating existing prognostic indices. ## MIPI
classic

```{r}
set.seed(102)
mcl_mipi_split <- initial_split(MCL, strata = status)
mcl_mipi_train <- training(mcl_mipi_split)
mcl_mipi_test <- testing(mcl_mipi_split)
mcl_mipi_folds <- vfold_cv(mcl_mipi_train, strata=status)
```

#Impute median values of MIPI

```{r}
mcl_mipi_recipe <- recipe(status~ calc_MIPI, data=mcl_mipi_train) %>% 
  step_impute_median(all_numeric_predictors()) 
 

mcl_mipi_recipe_prep <- prep(mcl_mipi_recipe)
juice(mcl_mipi_recipe_prep)
```

```{r}
mcl_mipi_spec <- logistic_reg() %>% 
  set_engine("glm")

mcl_mipi_spec
```

```{r}
mcl_mipi_wf <- workflow() %>% 
  add_recipe(mcl_mipi_recipe)

mcl_mipi_wf
```

```{r}
mcl_mipi_metrics <- metric_set(roc_auc, accuracy)

mcl_mipi_rs <- mcl_mipi_wf %>% 
  add_model(mcl_mipi_spec) %>% 
  fit_resamples(
    resamples=mcl_mipi_folds,
    metrics = mcl_mipi_metrics,
    control=control_resamples(save_pred = TRUE)
  )

mcl_mipi_rs
```

```{r}
collect_metrics(mcl_mipi_rs)
```

```{r}
mcl_mipi_final <- mcl_mipi_wf %>% 
  add_model(mcl_mipi_spec) %>% 
  last_fit(mcl_mipi_split)

mcl_mipi_final
collect_metrics(mcl_mipi_final)
```
```{r}
collect_predictions(mcl_mipi_final) %>% 
  conf_mat(status, .pred_class) 
```

#Univariate models incorporating existing prognostic indices. 
##Biological MIPI

```{r}
set.seed(102)
mcl_mipib_split <- initial_split(MCL, strata = status)
mcl_mipib_train <- training(mcl_mipib_split)
mcl_mipib_test <- testing(mcl_mipib_split)
mcl_mipib_folds <- vfold_cv(mcl_mipib_train, strata=status)
```

#Impute median values of MIPI

```{r}
mcl_mipib_recipe <- recipe(status~ biol_MIPI, data=mcl_mipib_train) %>% 
  step_impute_median(all_numeric_predictors()) 
 

mcl_mipib_recipe_prep <- prep(mcl_mipib_recipe)
juice(mcl_mipib_recipe_prep)
```

```{r}
mcl_mipib_spec <- logistic_reg() %>% 
  set_engine("glm")

mcl_mipib_spec
```

```{r}
mcl_mipib_wf <- workflow() %>% 
  add_recipe(mcl_mipib_recipe)

mcl_mipib_wf
```

```{r}
mcl_mipib_metrics <- metric_set(roc_auc, accuracy)

mcl_mipib_rs <- mcl_mipib_wf %>% 
  add_model(mcl_mipib_spec) %>% 
  fit_resamples(
    resamples=mcl_mipib_folds,
    metrics = mcl_mipib_metrics,
    control=control_resamples(save_pred = TRUE)
  )

mcl_mipib_rs
```


```{r}
collect_metrics(mcl_mipib_rs)
```

```{r}
mcl_mipib_final <- mcl_mipib_wf %>% 
  add_model(mcl_mipib_spec) %>% 
  last_fit(mcl_mipib_split)

mcl_mipib_final
collect_metrics(mcl_mipib_final)
```


```{r}
collect_predictions(mcl_mipib_final) %>% 
  conf_mat(status, .pred_class) 
```