COVID_endurance_impacts.qmd

---
title: "R Notebook"
author: "Félix Boudry"
date: "`r Sys.Date()`"
output:
    html_notebook:
        code_folding: hide
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  error = TRUE,
  message = FALSE,
  warning = FALSE,
  fig.width = 12,
  fig.height = 6
)

# Choosing seed for replicability
set.seed(123456)

# Import packages
library(rmarkdown)
library(tidyverse)
library(plyr)
library(janitor)
library(psych)
library(DTWBI)
library(stringi)
library(TraMineR)
library(reshape2)
library(cluster)
library(gtsummary)
library(ggpubr)
library(DT)
source(file = "Functions.R") # Import own functions
```

```{r Import}
# Import data sets
my_var <- "endurance_sport"
col_questions <- read.csv(file = "Data/col_question.csv")
analysis_data <- read.csv(file = "Data/Answers_completed.csv") %>% 
  `colnames<-`(col_questions$col_name) %>% 
  clean_names() %>% 
  remove_empty(which = c("rows", "cols")) %>% 
  select(-any_of("other_medication"))
```

```{r Compute}
# Compute new useful variables
analysis_data$duration_to_training_binary <-
  analysis_data$duration_to_training %>%
  stri_replace_all_regex(
    pattern = c(
      "Immediately|A few days later",
      ".*month.*|.*week.*"
    ),
    replacement = c("short", "long"),
    vectorize_all = FALSE
  )

analysis_data$difficulties_duration_binary <-
  analysis_data$difficulties_duration %>%
  stri_replace_all_regex(
    pattern = c("Nothing|1 week",
                ".*month.*|.*weeks.*"),
    replacement = c("short", "long"),
    vectorize_all = FALSE
  )

analysis_data$time_to_normal_training_volume_binary <-
  analysis_data$time_to_normal_training_volume %>%
  stri_replace_all_regex(
    pattern = c("1 week|1 month",
                ".*months.*|.*weeks.*"),
    replacement = c("short", "long"),
    vectorize_all = FALSE
  )

analysis_data["Training volume modifications"] <-
  analysis_data$modified_training_volume %>%
  stri_replace_all_regex(
    pattern = c(".*diminished.*",
                ".*augmented.*|.*No.*"),
    replacement = c("diminished", "augmented/unchanged"),
    vectorize_all = FALSE
  )

# Encode answers as labels
analysis_data_endcoded <-
  df_encode(input = analysis_data) # Create a data set with labeled data
```

```{r Statistics}
stat_results <- mapply(
  FUN = function(my_columns, my_colnames) {
    # Check number of possible answer to determine test
    answer_profile <-
      my_columns %>% unique() %>% na.omit() %>% length()
    usable_table <-
      na.omit(analysis_data[, c(my_var, my_colnames)])[[my_var]] %>%
      unique() %>% length()
    if (answer_profile < 2 && usable_table < 2) {
      return(NA)
    } else if (answer_profile >= 2 && usable_table > 1) {
      # Compute chi2 statistics
      chi2_stats(input = analysis_data,
                 var1 = my_colnames,
                 var2 = my_var)
    }
  },
  my_columns = analysis_data,
  my_colnames = names(analysis_data)
)
```

```{r Description}
# Analysis based on simple statistics and plots
descriptive_results <- mapply(
  FUN = function(my_columns, my_colnames) {
    # Check number of possible answer to determine plots
    answer_profile <-
      my_columns %>% unique() %>% na.omit() %>% length()
    usable_table <-
      na.omit(analysis_data[, c(my_var, my_colnames)])[[my_var]] %>%
      unique() %>% length()
    if (answer_profile < 2 && usable_table < 2) {
      # Only count answer if no groups
      my_count <-
        occurrences(input = analysis_data, feature = my_colnames)
      return(lst(my_count))
    } else if (answer_profile == 2 && usable_table > 1) {
      # Optimized plots for double bar plots
      my_count <-
        occurrences(input = analysis_data, feature = my_colnames)
      my_plot <-
        barplots(
          input = analysis_data,
          features = my_colnames,
          fill_feature = my_var,
          statistics_results = stat_results
        ) %>%
        getElement(name = my_colnames)
      return(lst(my_count, my_plot))
    } else if (answer_profile > 2 && usable_table > 1) {
      # Optimized plots for multiple bar plots
      my_count <-
        occurrences(input = analysis_data, feature = my_colnames)
      my_plot <-
        barplots(
          input = analysis_data,
          features = my_colnames,
          fill_feature = my_var,
          statistics_results = stat_results
        ) %>%
        getElement(name = my_colnames)
      return(lst(my_count, my_plot))
    } else {
      return(NA)
    }
  },
  my_columns = analysis_data,
  my_colnames = names(analysis_data)
)

descriptive_plots <- # Put all plots in a new list
  lapply(descriptive_results[!is.na(descriptive_results)], "[[", "my_plot") %>%
  discard(is.null)
descriptive_values <- # Put all counts in a new list
  lapply(descriptive_results[!is.na(descriptive_results)], "[[", "my_count") %>%
  discard(is.null)

age_18_25 <-
  round(count(analysis_data$age <= 25)[2, 2] / nrow(analysis_data) * 100,
        digits = 1)
age_26_59 <-
  round(
    count(26 <= analysis_data$age &
            analysis_data$age <= 59)[2, 2] / nrow(analysis_data) * 100,
    digits = 1
  )
age_60 <-
  round(count(analysis_data$age >= 60)[2, 2] / nrow(analysis_data) * 100,
        digits = 1)

pcr_test <-
  round(count(analysis_data$covid_test == "PCR test")[2, 2] / nrow(analysis_data) * 100,
        digits = 1)
antigenic_test <-
  round(
    count(analysis_data$covid_test == "Antigenic test")[2, 2] / nrow(analysis_data) * 100,
    digits = 1
  )

covid_type <- count(analysis_data$covid_type) %>% my_table()
multi_covid <-
  round(count(analysis_data$covid_multi)[2, 2] / nrow(analysis_data) * 100,
        digits = 1)

persistent_symptoms <-
  round(count(analysis_data$after_symptoms == "Yes")[2, 2] / nrow(analysis_data) *
          100,
        digits = 1)

after_fatigue_prev <-
  round(
    count(analysis_data$after_tired == "Yes")[2, 2] / count(analysis_data$after_symptoms == "Yes")[2, 2] *
      100,
    digits = 1
  )
after_breathlessness_prev <-
  round(
    count(analysis_data$after_breathlessness == "Yes")[2, 2] / count(analysis_data$after_symptoms == "Yes")[2, 2] *
      100,
    digits = 1
  )
after_cough_prev <-
  round(
    count(analysis_data$after_cough == "Yes")[2, 2] / count(analysis_data$after_symptoms == "Yes")[2, 2] *
      100,
    digits = 1
  )
after_headache_prev <-
  round(
    count(analysis_data$after_headache == "Yes")[2, 2] / count(analysis_data$after_symptoms == "Yes")[2, 2] *
      100,
    digits = 1
  )
after_muscle_prev <-
  round(
    count(analysis_data$after_muscular_pain == "Yes")[2, 2] / count(analysis_data$after_symptoms == "Yes")[2, 2] *
      100,
    digits = 1
  )

train_fatigue_prev <-
  round(count(analysis_data$training_tired == "Yes")[2, 2] / nrow(analysis_data) *
          100,
        digits = 1)
train_respi_prev <-
  round(
    count(analysis_data$training_respiratory_difficulties == "Yes")[2, 2] /
      nrow(analysis_data) * 100,
    digits = 1
  )
train_muscle_prev <-
  round(
    count(analysis_data$training_muscular_pain == "Yes")[2, 2] / nrow(analysis_data) *
      100,
    digits = 1
  )
train_concentration_prev <-
  round(
    count(analysis_data$training_concentration_difficulties == "Yes")[2, 2] /
      nrow(analysis_data) * 100,
    digits = 1
  )

complication_2w <-
  round(
    count(
      analysis_data$difficulties_duration == "2 weeks" |
        analysis_data$difficulties_duration == "1 week"
    )[2, 2] / nrow(analysis_data) * 100,
    digits = 1
  )
no_complications <-
  round(
    count(analysis_data$difficulties_duration == "Nothing")[2, 2] / nrow(analysis_data) *
      100,
    digits = 1
  )

data_end <- analysis_data[analysis_data$endurance_sport == "Yes", ]
data_nend <- analysis_data[analysis_data$endurance_sport == "No", ]

asthma_end <-
  round(nrow(data_end %>% filter(
    grepl("asthme|bie", pathology_description)
  )) / nrow(data_end) * 100, digits = 1)
asthme_nend <-
  round(nrow(data_nend %>% filter(
    grepl("asthme|bie", pathology_description)
  )) / nrow(data_nend) * 100, digits = 1)

data_symptoms <-
  analysis_data[, c(
    "endurance_sport",
    "tired",
    "headache",
    "cough",
    "runny_nose",
    "sore_throat",
    "fever",
    "muscular_pain"
  )] %>%
  `colnames<-`(
    c(
      "Endurance sport",
      "Unusual fatigue",
      "Headaches",
      "Cough",
      "Runny nose",
      "Sore throat",
      "Fever",
      "Muscle pain"
    )
  )

symptoms_plot <-
  data_symptoms %>%
  pivot_longer(cols = -"Endurance sport") %>%
  group_by("Endurance sport") %>%
  filter(value == "Yes") %>%
  ggplot(aes(fill = .$`Endurance sport`, x = name)) +
  geom_bar(position = "dodge") +
  geom_text(
    stat = "count",
    aes(label = ..count..),
    vjust = -0.2,
    position = position_dodge(width = 0.9)
  ) +
  labs(fill = "Practice an endurance sport") +
  ggtitle("Symptoms by population (END vs. NEND)") +
  scale_fill_manual(
    name = "",
    labels = c("NEND", "END"),
    values = c("darkgrey", "lightgrey")
  )
```

# Population description

## Anthropometrics

Presentation of the population studied.

In the population kept in this study there are
**`r descriptive_values$sex$woman`** womens and
**`r descriptive_values$sex$men`** mens and all of the
**`r nrow(analysis_data)`** subjects were infected by COVID-19.

```{r Anthropometrics}
analysis_data$covid_test %>%
  count() %>%
  my_table(col.names = c("Tests", "freq"))
describe(analysis_data[c("age", "height", "weight")], skew = FALSE) %>%
  round(digits = 1) %>% t() %>% my_table()
histograms(input = analysis_data,
           features = c("age", "height", "weight")) %>%
  ggarrange(plotlist = .,
            ncol = 3,
            nrow = 1)
```

## Training

```{r Training}
descriptive_plots[c("endurance_sport", "federal_license")] %>%
  ggarrange(plotlist = .,
            nrow = 1,
            ncol = 2)
walk(descriptive_plots[c("train_volume", "train_sessions_week", "train_method")], print)
count(df = analysis_data$what_sport)
```

In the studied population, **`r descriptive_values$endurance_sport$yes`** are endurance
trained athletes.

## Professions

Subject's professions are shown below:

```{r Professions}
count(df = analysis_data$profession)
```

## Pathologies

**`r count(analysis_data$pathology)[2, 2]`** have a known pathology.

```{r Pathologies}
count(df = analysis_data$pathology_description)
descriptive_plots$pathology_medication
```

# Answers by groups

## COVID

### Symptoms

```{r COVID symptoms}
descriptive_plots$symptoms_duration
symptoms_plot
descriptive_plots[c(
  "cough",
  "fever",
  "runny_nose",
  "headache",
  "sore_throat",
  "tired",
  "muscular_pain",
  "vomiting",
  "diarrhea",
  "breathlessness",
  "respiratory_difficulties",
  "other_symptoms"
)] %>%
  ggarrange(plotlist = .,
            ncol = 2,
            nrow = 2) %>%
  walk(print)
```

## Post COVID

### Symptoms

Binary : Up to a few days later = "short"; others =
"long"

```{r Post COVID symptoms}
descriptive_plots[c("after_symptoms", "duration_to_training_binary")] %>%
  ggarrange(plotlist = .,
            ncol = 2,
            nrow = 1)
descriptive_plots$duration_to_training
descriptive_plots[c(
  "after_cough",
  "after_runny_nose",
  "after_headache",
  "after_sore_throat",
  "after_tired",
  "after_muscular_pain",
  "after_vomiting",
  "after_diarrhea",
  "after_breathlessness",
  "after_respiratory_difficulties"
)] %>%
  ggarrange(plotlist = .,
            ncol = 2,
            nrow = 2) %>%
  walk(print)
```

### Impact during training

```{r Post COVID impact}
descriptive_plots[c(
  "training_tired",
  "training_respiratory_difficulties",
  "training_muscular_pain",
  "training_concentration_difficulties"
)] %>%
  ggarrange(plotlist = .,
            ncol = 2,
            nrow = 2)
descriptive_plots[c(
  "training_tired_note",
  "training_respiratory_difficulties_note",
  "training_muscular_pain_note",
  "training_concentration_difficulties_note"
)] %>% walk(print)
descriptive_plots$difficulties_duration_binary
descriptive_plots$difficulties_duration
descriptive_plots[c(
  "restart_training_medication",
  "respiratory_medication",
  "cortisone_medication",
  "vitamin_medication",
  "medication_duration"
)] %>%
  ggarrange(plotlist = .,
            ncol = 2,
            nrow = 2) %>%
  walk(print)
```

Binary : Up to 1 week = "short", others = "long"

### Training difficulties

```{r Training difficulties}
descriptive_plots[c(
  "high_intensity_training",
  "force_training",
  "endurance_training",
  "intermittent_training",
  "other_training"
)] %>%
  ggarrange(plotlist = .,
            ncol = 2,
            nrow = 2) %>%
  walk(print)
descriptive_plots[c(
  "longer_recuperation",
  "training_volume_back_normal",
  "modified_training_volume_binary"
)] %>%
  ggarrange(plotlist = .,
            ncol = 2,
            nrow = 2)
descriptive_plots[c(
  "modified_training_volume",
  "time_to_normal_training_volume",
  "endurance_training_more_difficult"
)] %>% walk(print)
```

Binary : Up to 1 month = "short"; others = "long"

### Hypoxia training

```{r Hypoxia training}
descriptive_plots[c(
  "hypoxia_training",
  "after_hypoxia_training",
  "x500_2000_hypoxia",
  "x2000_3000_hypoxia",
  "x3000_5500_hypoxia",
  "hypoxia_difficulties",
  "hypoxia_respiratory_difficulties"
)] %>%
  ggarrange(plotlist = .,
            ncol = 2,
            nrow = 2) %>%
  walk(print)
```

# Answer profiles

```{r Answer profiles}
my_profiles <- lst()
data_profiles <-
  analysis_data_endcoded$encoded_data %>%
  select_if(~ !any(is.na(.))) %>%
  select(# Removing columns with open questions and values
    -c(
      adult,
      profession,
      age,
      height,
      weight,
      pathology,
      pathology_description,
      pathology_medication,
      region,
      what_sport,
      train_volume,
      train_method,
      train_sessions_week,
      covid_positive,
      covid_test,
      covid_date,
      symptoms_duration,
      covid_type,
      covid_how_much,
      after_fever,
      medication_duration,
      restart_training,
      hypoxia_duration,
      difficulties_duration,
      duration_to_training,
      time_to_normal_training_volume
    )
  )
data_profiles_derived <- local.derivative.ddtw(data_profiles)

chains <-
  unite(data = data_profiles, col = "chains", sep = "-") %>%
  seqdef()
```

```{r Answer profiles plots}
matplot(t(data_profiles), type = 'l')
matplot(t(data_profiles_derived), type = 'l')

data_profiles %>%
  rownames_to_column() %>%
  gather(colname, value, -rowname) %>%
  ggplot() +
  geom_tile(aes(x = colname, y = rowname, fill = factor(value))) +
  scale_fill_brewer(palette = 10)

seqIplot(chains)

seqtrate(seqdata = chains, sel.states = NULL) %>%
  round(digits = 2) %>%
  melt() %>%
  ggplot(aes(Var2, Var1)) +
  geom_tile(aes(fill = value)) +
  geom_text(aes(label = round(value, 2))) +
  scale_fill_continuous(high = "#132B43",
                        low = "#56B1F7",
                        name = "Transitions")

seqsubm(chains, method = "TRATE") %>%
  round(digits = 2) %>%
  melt() %>%
  ggplot(aes(Var2, Var1)) +
  geom_tile(aes(fill = value)) +
  geom_text(aes(label = round(value, 2))) +
  scale_fill_continuous(high = "#132B43",
                        low = "#56B1F7",
                        name = "Substitution rate")

seqdist(
  chains,
  method = "OM",
  sm = seqsubm(chains, method = "TRATE"),
  indel = 1
) %>%
  round(digits = 2) %>%
  melt() %>%
  ggplot(aes(Var2, Var1)) +
  geom_tile(aes(fill = value)) +
  # geom_text(aes(label = round(value, 2))) +
  scale_fill_continuous(high = "#132B43",
                        low = "#56B1F7",
                        name = "Dissimilarities")

seqdist(
  chains,
  method = "OM",
  sm = seqsubm(chains, method = "TRATE"),
  indel = 1
) %>%
  agnes(diss = TRUE, method = "ward") %>%
  as.dendrogram() %>%
  plot(type = "rectangle")
```

# Unsupervised learning

```{r Unsupervised learning}
analysis_data$clust <- seqdist(
  chains,
  method = "OM",
  sm = seqsubm(chains, method = "TRATE"),
  indel = 1
) %>%
  agnes(diss = TRUE, method = "ward") %>%
  cutree(k = 2)
clust_summary <- tbl_summary(data = analysis_data, by = clust) %>%
  add_p() %>% as.data.frame()
datatable(clust_summary)
```

# Annexes

```{r Annexes}
questions <-
  lapply(col_questions,
         gsub,
         pattern = "_",
         replacement = " ") %>% as.data.frame()
datatable(questions)
```