associations.qmd

---
title: "Association Testing"
subtitle: "Comparing demographic data, clinical data, and outcomes across clusters."
author:
  - name: "Nathan Contantine-Cooke" 
    url: https://scholar.google.com/citations?user=2emHWR0AAAAJ&hl=en&oi=ao
    corresponding: true
    affiliations:
      - ref: HGU
      - ref: CGEM
  - name: "Karla Monterrubio-Gómez"
    url: https://scholar.google.com/citations?user=YmyxSXAAAAAJ&hl=en
    affiliations:
      - ref: HGU
  - name: "Riccardo E. Marioni"
    url: https://scholar.google.com/citations?hl=en&user=gA3Ik3MAAAAJ
    affiliations: 
      - ref: CGEM
  - name: "Catalina A. Vallejos"
    url: https://scholar.google.com/citations?user=lkdrwm0AAAAJ&hl=en&oi=ao
    affiliations:
      - ref: HGU
      - ref: Turing
#comments:
#  giscus: 
#    repo: quarto-dev/quarto-docs
---

## Introduction

```{R}
#| message: false
set.seed(123)
library(nnet) # Multinomial logistic regression
library(tidymodels) # Modelling framework
library(ranger) # Random forest
library(survival) # Kaplan-Meier
library(survminer) # Plot Kaplan-Meier curves
library(patchwork)
library(vip)
library(KONPsurv)

if (!require(plyr)) {
  install.packages("plyr")
}

########################
#-- Custom Functions --#
########################

#' @title Calculate descriptive statistics for non-normal continuous variables
#' @param vars Character vector of continuous variable names
#' @param data Dataframe object
#' @description
#' Finds the median, first quartile, third quartile, minimum and maximum
#'   statistics for non-normal continuous variables
#' @return A dataframe object with length(vars) rows and 7 columns
#' corresponding to the descriptive statistics calculated
describe_cont <- function(vars, data){
  if (tibble::is_tibble(data)) {
    data <- as.data.frame(data)
  }
  out_desc <- data.frame(Variable = character(),
                         n = double(),
                         Median = double(),
                         Q1 = double(),
                         Q3 = double(),
                         Min = double(),
                         Max = double()
  )
  for (i in 1:length(vars)) {
    var <- vars[i]
    med <- summary(data[, var])[c(3, 2, 5, 1, 6)]
    n <- sum(!is.na(data[, var]))
    out_desc[i, ] <- c(var, n, signif(med, 3))
  }
  out_desc
}
#' @title Calculate descriptive statistics for categorical variables
#' @param vars Character vector of categorical variable names
#' @param data Dataframe object
#' @description
#' Prints frequency and proportional tables for given categorical variables
#' @return A dataframe object with length(vars) rows and 3 columns
#' corresponding to the descriptive statistics calculated
describe_cat <- function(vars, data){
  if (tibble::is_tibble(data)) {
    data <- as.data.frame(data)
  }
  for (i in 1:length(vars)) {
    var <- vars[i]
    print(table(data[, var], useNA = "always"))
    print(prop.table(table(data[, var], useNA = "always")))
  }
}

kaplan_plot <- function(model, title, subtitle = TRUE, legend = TRUE){
  if (legend) {
    p <- ggsurvplot(model,
                    conf.int = TRUE,
                    surv.median.line = "none",
                    data = demographic,
                    pval = TRUE,
                    pval.method = TRUE,
                    ggtheme = theme_minimal() +
                      theme(plot.title = element_text(face = "bold",
                                                      size = 20)
                            ),
                    pval.size = 4,
                    legend = "bottom")
  } else {
        p <- ggsurvplot(model,
                    conf.int = TRUE,
                    surv.median.line = "none",
                    data = demographic,
                    pval = TRUE,
                    pval.method = TRUE,
                    ggtheme = theme_minimal() + 
                      theme(plot.title = element_text(face = "bold",
                                                      size = 20)
                            ),
                    pval.size = 4,
                    legend = "none")
  }
  if (subtitle) {
    p <- p + ggtitle(title, "Stratified by cluster") + xlab("Time (years)") 
  } else {
    p <- p + ggtitle(title) + xlab("Time (years)")
  }
  return(p)
} 
```

After [deeming the four-cluster LCMM to be the most suitable](selection.qmd), we
will now explore potential associations with cluster membership. We will
consider variables typically available at diagnosis, outcomes usually indicative
of a poor disease course, and treatments prescribed within one year of
diagnosis. For all categorical variables, frequency tables for the study
population and clusters have been generated in addition to either chi-squared
test or (if a cell has < 5 observations) Fisher's exact test results. For
continuous variables, the median, first quartile, and third quartile have been
reported for the study population and clusters in addition to ANOVA results. 

To streamline the analysis, descriptive functions have been created,
`describe_cat()` and `describe_cont()` which generate the aforementioned results

We also attempt to predict cluster membership based on these variables to
determine if cluster membership can be accurately predicted ahead of time using
these variables. We consider both multinomial logistic regression and random
forest predictive models for this purpose. 

## Characteristics available at diagnosis

We have found two variables available at diagnosis to be associated with cluster
membership: smoking at diagnosis ($p = 0.015$) and upper gastrointestinal
inflammation (Montreal L4; $p < 0.001$). 

```{R}
FCcumulative <- readRDS(paste0("/Volumes/igmm/cvallejo-predicct/cdi/processed/",
                               "FCcumulativeLongInc.RDS"))
model <- readRDS("cache/cubicbf.fits.RDS")[[4]]
model$pprob$class <- plyr::mapvalues(model$pprob$class,
                                     from = c(4, 3, 2, 1),
                                     to = c(2, 1, 3, 4))
demographic <- read.csv(paste0("/Volumes/igmm/cvallejo-predicct/cdi/data/",
                               "20220328/",
                               "Nathan_FCprogressioncohort.csv")
                         )
demographic <- subset(demographic, Number %in% model$pprob$id)
colnames(demographic)[1] <- "id"
classes <- model$pprob[, 1:2]
demographic <- merge(demographic, classes, by = "id")
demographic$cluster <- demographic$class
```

::: {.panel-tabset}

### Sex

#### Population 

```{R}
describe_cat("SEX", demographic)
```

#### Cluster = 1

```{R}
describe_cat("SEX", subset(demographic, class == 1))
```

#### Cluster = 2

```{R}
describe_cat("SEX", subset(demographic, class == 2))
```

#### Cluster = 3

```{R}
describe_cat("SEX", subset(demographic, class == 3))
```

#### Cluster = 4

```{R}
describe_cat("SEX", subset(demographic, class == 4))
```

#### Chi-squared test

```{R}
chisq.test(demographic$SEX, demographic$cluster)
```

### Age at diagnosis

#### Population 

```{R}
describe_cont("age", demographic)
```

#### Cluster = 1

```{R}
describe_cont("age", subset(demographic, class == 1))
```

#### Cluster = 2

```{R}
describe_cont("age", subset(demographic, class == 2))
```

#### Cluster = 3

```{R}
describe_cont("age", subset(demographic, class == 3))
```

#### Cluster = 4

```{R}
describe_cont("age", subset(demographic, class == 4))
```

#### ANOVA

```{R}
summary(aov(age ~ class, data = demographic))
```

### Smoking at diagnosis

If a participant has reported to have not smoked in the past, then it is assumed
the participant did not smoke at IBD diagnosis. Otherwise smoking status is
taken from the participant's response to being asked if they were smoking when
diagnosed with IBD. 

```{R}
colnames(demographic)[colnames(demographic) == "SMOKER..Y.N."] <- "smoking"
demographic$smoking <- factor(demographic$smoking, labels = c("no", "yes"))
```

#### Population

```{R}
describe_cat("smoking", demographic)
```

#### Cluster = 1 

```{R}
describe_cat("smoking", subset(demographic, class == 1))
```

#### Cluster = 2

```{R}
describe_cat("smoking", subset(demographic, class == 2))
```

#### Cluster = 3

```{R}
describe_cat("smoking", subset(demographic, class == 3))
```

#### Cluster = 4

```{R}
describe_cat("smoking", subset(demographic, class == 4))
```


#### Fisher's exact test

```{R}
fisher.test(demographic$cluster, demographic$smoking)
```

### Diagnostic FCAL

#### Population 

```{R}
FCcumulative <- readRDS(paste0("/Volumes/igmm/cvallejo-predicct/processed-data/cdi/",
                               "FCcumulativeLongInc.RDS"))
FCcumulative <- FCcumulative[order(FCcumulative$time), ]
DiagFC <- rep(NA, nrow(demographic))
for (i in 1:nrow(demographic)) {
  subject <- demographic[i, "id"]
  subject.data <- subset(FCcumulative, id == subject)
  DiagFC[i] <- subject.data[1, "value"] # Already sorted by time
}
demographic$DiagFC <- DiagFC

describe_cont("DiagFC", demographic)
```

#### Cluster = 1

```{R}
describe_cont("DiagFC", subset(demographic, class == 1))
```

#### Cluster = 2

```{R}
describe_cont("DiagFC", subset(demographic, class == 2))
```

#### Cluster = 3

```{R}
describe_cont("DiagFC", subset(demographic, class == 3))
```

#### Cluster = 4

```{R}
describe_cont("DiagFC", subset(demographic, class == 4))
```

#### ANOVA

```{R}
summary(aov(DiagFC ~ class, data = demographic))
```

### Montreal behaviour 

```{R}
colnames(demographic)[colnames(demographic) == "BEHAVIOUR.AT.DIAGNOSIS"] <- "behaviour"
demographic$behaviour <- as.factor(demographic$behaviour)
```

#### Population

```{R}
describe_cat("behaviour", demographic)
```

#### Cluster = 1 

```{R}
describe_cat("behaviour", subset(demographic, class == 1))
```

#### Cluster = 2

```{R}
describe_cat("behaviour", subset(demographic, class == 2))
```

#### Cluster = 3

```{R}
describe_cat("behaviour", subset(demographic, class == 3))
```

#### Cluster = 4

```{R}
describe_cat("behaviour", subset(demographic, class == 4))
```

#### Fisher's exact test

```{R}
fisher.test(demographic$cluster, demographic$behaviour)
```

### Perianal

```{R}
colnames(demographic)[colnames(demographic) == "PERIANAL.DISEASE..Y.N."] <- "peri"
demographic$peri <- factor(demographic$peri, labels = c("no", "yes"))
```

#### Population

```{R}
describe_cat("peri", demographic)
```

#### Cluster = 1 

```{R}
describe_cat("peri", subset(demographic, class == 1))
```

#### Cluster = 2

```{R}
describe_cat("peri", subset(demographic, class == 2))
```

#### Cluster = 3

```{R}
describe_cat("peri", subset(demographic, class == 3))
```

#### Cluster = 4

```{R}
describe_cat("peri", subset(demographic, class == 4))
```

#### Fisher's exact test

```{R}
fisher.test(demographic$cluster, demographic$peri)
```

### Montreal location

```{R}
demographic$LOCATION <- factor(demographic$LOCATION, labels = c("L1", "L2", "L3"))
```

#### Population

```{R}
describe_cat("LOCATION", demographic)
```

#### Cluster = 1 

```{R}
describe_cat("LOCATION", subset(demographic, class == 1))
```


#### Cluster = 2

```{R}
describe_cat("LOCATION", subset(demographic, class == 2))
```

#### Cluster = 3 

```{R}
describe_cat("LOCATION", subset(demographic, class == 3))
```

##### Cluster = 4

```{R}
describe_cat("LOCATION", subset(demographic, class == 4))
```


#### Fisher's exact test

```{R}
fisher.test(demographic$cluster, demographic$LOCATION, workspace = 800000)
```

### Upper GI

```{R}
colnames(demographic)[colnames(demographic) == "L4.Modifier"] <- "L4"
demographic$L4 <- factor(demographic$L4, labels = c("no", "yes"))
```

#### Population

```{R}
describe_cat("L4", demographic)
```

#### Cluster = 1 

```{R}
describe_cat("L4", subset(demographic, class == 1))
```

#### Cluster = 2

```{R}
describe_cat("L4", subset(demographic, class == 2))
```

#### Cluster = 3

```{R}
describe_cat("L4", subset(demographic, class == 3))
```

#### Cluster = 4

```{R}
describe_cat("L4", subset(demographic, class == 4))
```

#### Fisher's exact test

```{R}
fisher.test(demographic$cluster, demographic$L4)
```

:::

## Associations with treatment information

To explore the potential presence of treatment effects, we have tested for
associations between cluster membership and common IBD prescriptions. All 
treatments were prescribed within one year of diagnosis unless otherwise stated.
Thioprine use ($p = 0.02$), biologic use within three months of diagnosis
($p < 0.001$), and biologic use within one year of diagnosis ($p < 0.001$) were
found to be significantly associated with cluster membership. 

::: {.panel-tabset}

### 5ASA

#### Population

```{R}
describe_cat("X5ASA", demographic)
```

#### Cluster = 1

```{R}
describe_cat("X5ASA", subset(demographic, class == 1))
```

#### Cluster = 2

```{R}
describe_cat("X5ASA", subset(demographic, class == 2))
```

#### Cluster = 3

```{R}
describe_cat("X5ASA", subset(demographic, class == 3))
```

#### Cluster = 4

```{R}
describe_cat("X5ASA", subset(demographic, class == 4))
```

#### Fisher's exact test

```{R}
fisher.test(demographic$cluster, demographic$X5ASA)
```

### Thioprine

```{R}
colnames(demographic)[colnames(demographic) == "THIOPURINE.WITHIN.1st.YEAR..Y.N."] <- "use_thio"
```

#### Population

```{R}
describe_cat("use_thio", demographic)
```

#### Cluster = 1

```{R}
describe_cat("use_thio", subset(demographic, class == 1))
```

#### Cluster = 2

```{R}
describe_cat("use_thio", subset(demographic, class == 2))
```

#### Cluster = 3

```{R}
describe_cat("use_thio", subset(demographic, class == 3))
```

#### Cluster = 4

```{R}
describe_cat("use_thio", subset(demographic, class == 4))
```

#### Fisher's exact test

```{R}
fisher.test(demographic$cluster, demographic$use_thio)
```

### Corticosteroids

```{R}
colnames(demographic)[colnames(demographic) == "CORTICOSTEROIDS.WITHIN.1st.YEAR...Y.N."] <- "use_cortico"
demographic$use_cortico <- plyr::mapvalues(demographic$use_cortico, "y", "Y")
```

#### Population

```{R}
describe_cat("use_cortico", demographic)
```

#### Cluster = 1

```{R}
describe_cat("use_cortico", subset(demographic, class == 1))
```

#### Cluster = 2

```{R}
describe_cat("use_cortico", subset(demographic, class == 2))
```

#### Cluster = 3

```{R}
describe_cat("use_cortico", subset(demographic, class == 3))
```

#### Cluster = 4

```{R}
describe_cat("use_cortico", subset(demographic, class == 4))
```

#### Fisher's exact test

```{R}
fisher.test(demographic$cluster, demographic$use_cortico)
```

### Methotrexate

```{R}
colnames(demographic)[colnames(demographic) == "MTX.WITHIN.1st.YEAR..Y.N."] <- "use_mtx"
```

#### Population

```{R}
describe_cat("use_mtx", demographic)
```

#### Cluster = 1

```{R}
describe_cat("use_mtx", subset(demographic, class == 1))
```

#### Cluster = 2

```{R}
describe_cat("use_mtx", subset(demographic, class == 2))
```

#### Cluster = 3

```{R}
describe_cat("use_mtx", subset(demographic, class == 3))
```

#### Cluster = 4

```{R}
describe_cat("use_mtx", subset(demographic, class == 4))
```

#### Fisher's exact test

```{R}
fisher.test(demographic$cluster, demographic$use_mtx)
```

### EEN

```{R}
colnames(demographic)[colnames(demographic) == "EEN..Y.N."] <- "EEN"
demographic$EEN <- na_if(demographic$EEN, "")
```

#### Population

```{R}
describe_cat("EEN", demographic)
```

#### Cluster = 1

```{R}
describe_cat("EEN", subset(demographic, class == 1))
```

#### Cluster = 2

```{R}
describe_cat("EEN", subset(demographic, class == 2))
```

#### Cluster = 3

```{R}
describe_cat("EEN", subset(demographic, class == 3))
```

#### Cluster = 4

```{R}
describe_cat("EEN", subset(demographic, class == 4))
```

#### Fisher's exact test

```{R}
fisher.test(demographic$cluster, demographic$EEN)
```

### Biologic use within three months

```{R}
colnames(demographic)[colnames(demographic) == "BIOLOGIC.MONO...COMBO.WITHIN.3M"] <- "use_bio3"
```

#### Population

```{R}
describe_cat("use_bio3", demographic)
```

#### Cluster = 1

```{R}
describe_cat("use_bio3", subset(demographic, class == 1))
```

#### Cluster = 2

```{R}
describe_cat("use_bio3", subset(demographic, class == 2))
```

#### Cluster = 3

```{R}
describe_cat("use_bio3", subset(demographic, class == 3))
```

#### Cluster = 4

```{R}
describe_cat("use_bio3", subset(demographic, class == 4))
```

#### Fisher's exact test

```{R}
fisher.test(demographic$cluster, demographic$use_bio3)
```

### Biologic use within one year

```{R}
colnames(demographic)[colnames(demographic) == "BIOLOGIC.WITHIN.1st.YEAR..Y.N."] <- "use_bio"
```

#### Population

```{R}
describe_cat("use_bio", demographic)
```

#### Cluster = 1

```{R}
describe_cat("use_bio", subset(demographic, class == 1))
```

#### Cluster = 2

```{R}
describe_cat("use_bio", subset(demographic, class == 2))
```

#### Cluster = 3

```{R}
describe_cat("use_bio", subset(demographic, class == 3))
```

#### Cluster = 4

```{R}
describe_cat("use_bio", subset(demographic, class == 4))
```

#### Fisher's exact test

```{R}
fisher.test(demographic$cluster, demographic$use_bio)
```

### Date of diagnosis

Here, we consider the number of days from the start of the study period
(2005-01-1) to when a subject was diagnosed to determine if when a subject was
diagnosed has a statistically significant impact on which cluster a subject was
assigned to. 

```{R}
demographic <- datefixR::fix_dates(demographic, "DATE.OF.DIAGNOSIS")
demographic$Days.From.Start <- as.numeric(demographic$DATE.OF.DIAGNOSIS) -
  as.numeric(as.Date("2005-01-01"))

describe_cont("Days.From.Start", demographic)
```

#### Cluster = 1

```{R}
describe_cont("Days.From.Start", subset(demographic, class == 1))
```

#### Cluster = 2

```{R}
describe_cont("Days.From.Start", subset(demographic, class == 2))
```

#### Cluster = 3

```{R}
describe_cont("Days.From.Start", subset(demographic, class == 3))
```

#### Cluster = 4

```{R}
describe_cont("Days.From.Start", subset(demographic, class == 4))
```

#### ANOVA

```{R}
summary(aov(Days.From.Start ~ class, data = demographic))
```


:::

## Predictive models

Here, we consider whether we can accurately predict cluster membership using age,
sex, and the variables we found to be significantly associated with
cluster membership: smoking at diagnosis, upper gastrointestinal inflammation
(Montreal L4), and being prescribed a biologic therapeutic within one year of
diagnosis. We consider two types of models: multinomial logistic
regression models and a random forest classifiers using the `{nnet}` and
`{ranger}` packages. We use the `{tidymodels}` framework for fitting these
models.

The data are split into a 75:25 train:test split and k-fold cross validation 
is used. Reported metrics and the confusion matrix are based on the test data
split. 


We attempt to predict cluster membership for all four clusters.

```{R}
demographic$cluster <- as.factor(demographic$cluster)
demographic$cluster = relevel(demographic$cluster, ref = 2) # 

data_split <- initial_split(demographic, prop = 0.75, strata = cluster)
train_data <- training(data_split)
test_data  <- testing(data_split)

folds <- vfold_cv(train_data, v = 4, strata = cluster)

class_rec <- recipe(cluster ~ age +
                      SEX  +
                      smoking +
                      LOCATION +
                      L4 +
                      behaviour +
                      peri +
                      DiagFC +
                      use_bio +
                      use_bio3,
                    data = test_data)
```

### Multinomial logistic regression

```{R}
mlr_mod <- multinom_reg(penalty = tune(),
                        mode = "classification") %>%
  set_engine("nnet")

class_wflow <- 
  workflow() %>% 
  add_model(mlr_mod) %>% 
  add_recipe(class_rec)

mlr_grid <- grid_regular(
                         penalty(),
                         levels = 5)

class_fit <- class_wflow %>% 
    tune_grid(
    resamples = folds,
    grid = mlr_grid
    )

best_mlr <- class_fit %>%
  select_best("accuracy")

final_wf <- 
  class_wflow %>% 
  finalize_workflow(best_mlr)


final_fit <- 
  final_wf %>%
  last_fit(data_split) 

temp <- final_fit %>% extract_fit_engine()

temp$call <- quote(nnet::multinom(formula = cluster ~ age +
                                      SEX  +
                                      smoking +
                                      LOCATION +
                                      L4 +
                                      behaviour +
                                      peri +
                                      DiagFC +
                                      use_bio +
                                      use_bio3,
                                  data = train_data,
                                  decay = ~0.00316227766016838, 
                                  trace = FALSE))

hmm <- summary(temp)
```

#### Odd ratios 

Point estimates:

```{R}
knitr::kable(round(exp(hmm$coefficients), 2)) 
```

Lower 95% confidence interval:

```{R}
knitr::kable(round(exp(hmm$coefficients - qnorm(0.975) * hmm$standard.errors), 2))
```

Upper 95% confidence interval:

```{R}
knitr::kable(round(exp(hmm$coefficients + qnorm(0.975) * hmm$standard.errors), 2))
```

#### Metrics

```{R}
knitr::kable(
  final_fit %>%
    collect_metrics(), 
  col.names = c("Metric", "Estimator", "Estimate", "Config"),
  align = "cccc") 
```

#### Confusion matrix

```{R}
results <- final_fit %>%
  collect_predictions() 

knitr::kable(table(results$cluster, results$.pred_class), row.names = TRUE)
```

### Random forest

#### Metrics

```{R}
rf_mod <- rand_forest(mtry = tune(), 
                      trees = tune(),
                      mode = "classification") %>%
  set_engine("ranger",
             num.threads = parallel::detectCores(),
             importance = "permutation")

class_wflow <- 
  workflow() %>% 
  add_model(rf_mod) %>% 
  add_recipe(class_rec)

tree_grid <- grid_regular(mtry(range(c(1,8))),
                          trees(),
                          levels = 5)

class_fit <- class_wflow %>% 
    tune_grid(
    resamples = folds,
    grid = tree_grid
    )

best_tree <- class_fit %>%
  select_best("accuracy")

final_wf <- 
  class_wflow %>% 
  finalize_workflow(best_tree)

final_fit <- 
  final_wf %>%
  last_fit(data_split) 

knitr::kable(
  final_fit %>% 
    extract_fit_parsnip() %>% 
    vi()
)

knitr::kable(
  final_fit %>%
    collect_metrics(), 
  col.names = c("Metric", "Estimator", "Estimate", "Config"),
  align = "cccc") 
```

#### Confusion matrix

```{R}
results <- final_fit %>%
  collect_predictions() 
knitr::kable(table(results$cluster, results$.pred_class), row.names = TRUE)
```

Model performance for both the multinomial logistic regression model and the
random forest classifier is very poor. 

## Session information

```{R}
#| echo: false
pander::pander(sessionInfo())
```

## {.appendix}

<img class ="center" src="images/MRC_HGU_Edinburgh RGB.png" alt="MRC Human Genetics Unit logo" height = 50px>
<img src="images/cgem-logo.png" alt="Centre for Genomic & Experimental Medicine logo" height = 50px> 


## Acknowledgments {.appendix}

This work is funded by the Medical Research Council & University of Edinburgh
via a Precision Medicine PhD studentship (MR/N013166/1, to **NC-C**) 

## Author contributions {.appendix} 

**NC-C** wrote the analysis. **KM** and **CAV** performed code review and
contributed suggestions. **RM** contributed functions for calculating
descriptive statistics. **KM**, **RM** and **CAV** provided feedback. 

## Reuse {.appendix}

Licensed by 
<a href="https://creativecommons.org/licenses/by/4.0/">CC BY</a>
except for the MRC Human Genetics Unit, The University of Edinburgh, and Centre for Genomic & Experimental Medicine logos or unless otherwise stated.