MWF_Models.Rmd

---
title: "MWF Models_PASE"
author: "Narlon Cassio"
date: "February 09, 2022"
output: 
  pdf_document: 
    toc: yes
    number_sections: yes
    toc_depth: 5
  word_document: 
    toc: yes
    toc_depth: 5
geometry: "left = 1cm, right = 1cm, top = 1cm, bottom = 2.5cm"
header-includes:
- \usepackage{caption}
- \usepackage{pdflscape}
- \newcommand{\blandscape}{\begin{landscape}}
- \newcommand{\elandscape}{\end{landscape}}
editor_options: 
  chunk_output_type: console
---

\newpage
```{r setup, include = FALSE}
knitr::opts_chunk$set(
  message = FALSE,
  fig.height = 5, fig.width = 6, fig.align="center"
)
```
\small
# Loading packages
```{r, }
#------------------------------------------------------------------#
# Loading packages                                              ####
#------------------------------------------------------------------#
library(tidyverse)
library(tableone)
library(gvlma)
library(openxlsx)
library(broom)
library(psych)
library(knitr)
library(rstatix)
library(performance)
```

# Loading data
```{r, message = FALSE}
#------------------------------------------------------------------#
# Loading data                                                  ####
#------------------------------------------------------------------#
# Outcome and descriptive data
all_data_clean <- read.xlsx("all_data_clean_final.xlsx")

```
# Data analysis
## Descriptives
```{r, }
# Demographics
kable(digits = 2, caption = "Demographics",
      describe(all_data_clean[c("age", 
                                "height", 
                                "weight",
                                "bmi",
                                "moca",
                                "mmse",
                                "overall_fall_risk_score",
                                "meters_walked",
                                "pase")], fast = TRUE, IQR = TRUE))


## Categorical demographics data 
kable(print(printToggle = FALSE, CreateTableOne(data = all_data_clean, c("sex", "fazekas_score", "education"))))
kable(print(printToggle = FALSE, CreateTableOne(data = all_data_clean, 
                                                c("fci_1_arthritis",
                                                  "fci_2_osteoporosis",
                                                  "fci_3_asthma",
                                                  "fci_4_copd_ards_or_emphysema",
                                                  "fci_5_angina",
                                                  "fci_6_congestive_heart_failure_or_heart_disease",
                                                  "fci_7_heart_attack_myocardial_infarct",
                                                  "fci_8_neurological_disease",
                                                  "fci_9_stroke_or_tia",
                                                  "fci_10_peripheral_vascular_disease",
                                                  "fci_11_diabetes_type_i_and_ii",
                                                  "fci_12_upper_gastrointestinal_disease",
                                                  "fci_13_depression",                                
                                                  "fci_14_anxiety_or_panic_disorders",
                                                  "fci_15_visual_impairment",
                                                  "fci_16_hearing_impairment",
                                                  "fci_17_degenerative_disc_disease",
                                                  "fci_18_obesity_and_or_body_mass_index_30",
                                                  "fci_19_thyroid_disease",
                                                  "fci_20_cancer",
                                                  "fci_21_hypertension",
                                                  "fci_total"))))
## Total Functional commorbidity index
kable(print(printToggle = FALSE, CreateTableOne(data = all_data_clean, c("fci_total"), 
                           includeNA = TRUE), nonnormal = c("fci_total")))

# Myelin water fraction
kable(digits = 2, caption = "Myelin water fraction",
      describe(all_data_clean[c("Ant_CR_mean",
                             "Post_CR_mean",
                             "Sup_CR_mean",
                             "BCC_mean",
                             "GCC_mean",
                             "SCC_mean",
                             "Cingulum_mean", 
                             "Ext_Cap_mean", 
                             "Fornix_mean",
                             "Sag_Strat_mean",
                             "Ant_IC_mean",
                             "Post_IC_mean", 
                             "Retro_IC_mean",
                             "Post_TR_mean",
                             "SLF_mean",
                             "MWF_TOT_WM_mean",
                             "MWF_TOT_WM_ero_mean",
                             "wmh_cm3",
                             "wmh_cm3_lg",
                             "eicv_cm3")], fast = TRUE, IQR = TRUE))

# Other structural data
kable(print(printToggle = FALSE, CreateTableOne(data = all_data_clean, c("wmh_cm3", "wmh_cm3_lg","eicv_cm3"), 
                           includeNA = TRUE), nonnormal = c("wmh_cm3")))

```
### Comparing participants between studies

```{r }
all_data_clean <- all_data_clean %>% 
  mutate(cohort = ifelse(str_detect(id, "FALLERS2_..."), "CogMob2", "RVCI")) %>% 
  mutate(scanner = ifelse(row_number() <= 72, "old","new"))

kable(print(printToggle = FALSE, CreateTableOne(
  data = all_data_clean, c(
    "age","height", "weight",
    "sex", "fazekas_score", "education",
    "bmi", "moca", "mmse","fci_total",
    "fci_1_arthritis",
    "fci_2_osteoporosis",
    "fci_3_asthma",
    "fci_4_copd_ards_or_emphysema",
    "fci_5_angina",
    "fci_6_congestive_heart_failure_or_heart_disease",
    "fci_7_heart_attack_myocardial_infarct",
    "fci_8_neurological_disease",
    "fci_9_stroke_or_tia",
    "fci_10_peripheral_vascular_disease",
    "fci_11_diabetes_type_i_and_ii",
    "fci_12_upper_gastrointestinal_disease",
    "fci_13_depression",                                
    "fci_14_anxiety_or_panic_disorders",
    "fci_15_visual_impairment",
    "fci_16_hearing_impairment",
    "fci_17_degenerative_disc_disease",
    "fci_18_obesity_and_or_body_mass_index_30",
    "fci_19_thyroid_disease",
    "fci_20_cancer",
    "fci_21_hypertension"), 
  includeNA = TRUE, test = TRUE, strata = "cohort"), explain = TRUE, noSpaces = TRUE,
  nonnormal = c("bmi", "moca", "mmse","fci_total")))

```

## Linear models
```{r, }
#------------------------------------------------------------------#
# Regression models                                             ####
#------------------------------------------------------------------#

# Function to summarize lm output
lm_results <- function(x) {
  my_list <- list()
  my_list$Summary <- (summary(x))
  my_list$Assumtions <- (gvlma(x))
  return(my_list)
}

# Creating contrasts for multi-level categorical variables
all_data_clean <- within(all_data_clean, {
  fazekas_c <- C(as.factor(fazekas_score), treatment, 2)
  education_c <- C(as.factor(education), treatment, 3)
  print(attributes(fazekas_c))
  print(attributes(education_c))
})

```
\newpage

## Models
### PASE
```{r, fig.height = 10, fig.width = 9, fig.align="center" }
## Select variable names for models
outcome_var <- colnames(all_data_clean[c(39:53,56,57)])
outcome_var %>% kable(col.names = "Outcome variables")

# List to store models
pase_lm = list()

# Run models in loop script
for (var in outcome_var){
  # Formula
  form <- paste(var, "~","pase")

    # Print results to the list
  pase_lm[var] <- list(lm(form,
          data = all_data_clean))
}

## Printing models
for (model in pase_lm){
  print(formula(model), showEnv = FALSE)
  print(lm_results(model))
}
```
### PASE + covariates
```{r, fig.height = 10, fig.width = 9, fig.align="center" }
## Select variable names for models
outcome_var <- colnames(all_data_clean[c(39:53,56,57)])
outcome_var %>% kable(col.names = "Outcome variables")

pase_adj_lm = list()

# Run models in loop script
for (var in outcome_var){
  # Formula
  form <- paste(var,
    "~","pase", 
    "+ age + eicv_cm3 + sex + bmi + education_c") # Covariates

    # Print results to the list
  pase_adj_lm[var] <- list(lm(form,
          data = all_data_clean))
}

## Printing models
for (model in pase_adj_lm){
  print(formula(model), showEnv = FALSE)
  print(lm_results(model))
  print(check_model(model))
}
```

### PASE + covariates + WMH
```{r, fig.height = 10, fig.width = 9, fig.align="center" }
pase_adj_wmh_lm <- list(lm(Ant_CR_mean ~ pase
                                + age
                                + eicv_cm3
                                + bmi
                                + sex
                                + education_c
                                + wmh_cm3_lg,
                                data = all_data_clean),
                    
                      lm(GCC_mean ~ pase
                                + age
                                + eicv_cm3
                                + bmi
                                + sex
                                + education_c
                                + wmh_cm3_lg,
                                data = all_data_clean),
                      
                      lm(Sag_Strat_mean ~ pase
                                + age
                                + eicv_cm3
                                + bmi
                                + sex
                                + education_c
                                + wmh_cm3_lg,
                                data = all_data_clean),
                     
                     lm(MWF_TOT_WM_ero_mean ~ pase
                                + age
                                + eicv_cm3
                                + bmi
                                + sex
                                + education_c
                                + wmh_cm3_lg,
                                data = all_data_clean),
                        
                        lm(MWF_TOT_WM_mean ~ pase
                                + age
                                + eicv_cm3
                                + bmi
                                + sex
                                + education_c
                                + wmh_cm3_lg,
                                data = all_data_clean))

for (model in pase_adj_wmh_lm){
  print(lm_results(model))
  print(check_model(model))
}
```
# Results
## Model summaries
```{r }
## Unadjusted
all_models_unadjusted_var = tibble()

for (model in pase_lm){
  output <- rbind(tidy(model, conf.int = TRUE, conf.level = 0.95))[2,]
  output$outcome <-  as.character(model$terms[[2]]) 
  all_models_unadjusted_var <- rbind(all_models_unadjusted_var,output)
}

## Printing estimates
all_models_unadjusted_var %>% 
  arrange(p.value) %>% 
  kable(digits = 3)

```

```{r}
## Adjusted
all_models_adjusted_var = tibble()

for (model in pase_adj_lm){
  output <- rbind(tidy(model, conf.int = TRUE, conf.level = 0.95))[2,]
  output$outcome <-  as.character(model$terms[[2]]) 
  all_models_adjusted_var <- rbind(all_models_adjusted_var,output)
}

## Printing estimates
all_models_adjusted_var %>% 
  arrange(p.value) %>% 
  kable(digits = 3)
```

```{r}
## Adjusted (WMH)
all_models_adjusted_var_wmh = tibble()

for (model in pase_adj_wmh_lm){
  output <- rbind(tidy(model, conf.int = TRUE, conf.level = 0.95))[2,]
  output$outcome <-  as.character(model$terms[[2]]) 
  all_models_adjusted_var_wmh <- rbind(all_models_adjusted_var_wmh,output)
}

## Printing estimates
all_models_adjusted_var_wmh %>% 
  arrange(p.value) %>% 
  kable(digits = 3)
```

## Plotting results
### PASE models
```{r, fig.height = 10, fig.width = 4.5, fig.align="center" }
## Prepting data for ploting (MWF by ROI)
all_data_clean %>% 
  pivot_longer(names_to = "roi", values_to = "mwf", cols = c(39:57)) %>% 
  select(id, roi, mwf, pase) %>% 
  mutate(roi = str_replace_all(roi, c("Ant_CR_mean" = "Corona Radiata, anterior",
                                      "GCC_mean" = "Corpus Callosum, genu",
                                      "Sag_Strat_mean" = "Sagittal Stratum"))) %>%
  filter(roi == "Corona Radiata, anterior" |
         roi == "Corpus Callosum, genu"|
         roi == "Sagittal Stratum") %>% 
mutate(anat_loc = ifelse(roi == "Corona Radiata, anterior", 3,
                  ifelse(roi == "Corpus Callosum, genu", 2,
                  ifelse(roi == "Sagittal Stratum", 1, NA)))) %>% 
  filter(is.na(pase)==FALSE) %>% 

## Plotting MWF ROI data
ggplot(aes(pase, mwf, colour = as.factor(anat_loc))) +
  geom_point(alpha = 0.8, size = 2.5) +
  scale_colour_manual(values = c("#A22917", # orange
                                 "#DF5E28", # light green
                                 "#EDB140")) +
  xlab("\n Physical activity levels") +
  ylab("Myelin water fraction (%) \n") +
  scale_y_continuous(breaks = seq(0.00, 30, by = 2.0)) +
  scale_x_continuous(breaks = seq(0, 300, by = 50)) +
  facet_wrap(~reorder(roi, anat_loc), ncol = 1) +
  geom_smooth(alpha = 0.3, method = lm, se = TRUE) +
  theme_minimal() +
  theme(legend.position = "none",
        strip.text = element_text(size = 12),
        axis.text = element_text(colour = "black", size = 12),
        axis.title = element_text(colour = "black", size = 12),
        panel.grid.major = element_blank())
ggsave("plots/mwf_pase.jpeg", height = 10, width = 4.5, dpi = 500)

```

```{r, fig.height = 4, fig.width = 5.5}
## Prepting data for ploting whole-brain white matter
all_data_clean %>% 
  pivot_longer(names_to = "roi", values_to = "mwf", cols = c(39:57)) %>% 
  select(id, roi, mwf, pase) %>% 
  mutate(roi = str_replace_all(roi, c("MWF_TOT_WM_ero_mean" = "Whole-brain white matter"))) %>%
  filter(roi == "Whole-brain white matter") %>%
  filter(is.na(pase)==FALSE) %>% 

## Plotting MWF whole-brain white matter data
ggplot(aes(pase, mwf)) +
  geom_point(alpha = 0.7, colour = "black",  size = 2.5) +
  #scale_colour_gradient(low = "#48caa2", high = "#f57100") +
  xlab("\n Physical activity levels") +
  ylab("Myelin water fraction (%) \n") +
  scale_y_continuous(breaks = seq(0.00, 30, by = 1), n.breaks = 10) +
  scale_x_continuous(breaks = seq(0, 300, by = 50)) +
  facet_wrap(~roi, ncol = 1) +
  geom_smooth(colour = "black", alpha = 0.3, method = lm, se = TRUE) +
  theme_minimal() +
  theme(legend.position = "none",
        strip.text = element_text(size = 8),
        axis.text = element_text(colour = "black", size = 8),
        axis.title = element_text(colour = "black", size = 8),
        panel.grid.major = element_blank())
ggsave("plots/mwf_whole-brain_pase.jpeg", height = 3, width = 4.5, dpi = 500)

all_data_clean %>% 
  ggplot(aes(MWF_TOT_WM_ero_mean, MWF_TOT_WM_mean)) +
  geom_point(size = 2, alpha = 0.75) +
  ylab("Whole-brain white matter, original") +
  xlab("Whole-brain white matter, eroded") +
  theme_minimal()
ggsave("plots/mwf_whole-brain_eroded_vs_original.jpeg", height = 3, width = 4)

```