sleep_thickness_models.Rmd

---
title: "Sleep, Cortical Thickness, and Cognition Analysis"
author: "Nárlon Cássio Boa Sorte Silva, PhD"
date: "July 07, 2022"
output: 
  pdf_document: 
    toc: yes
    number_sections: yes
    toc_depth: 6
  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.align = 'center')
```

\small

# Loading packages
```{r, }
#------------------------------------------------------------------#
# Loading packages                                              ####
#------------------------------------------------------------------#
library(tidyverse)
library(tableone)
library(gvlma)
library(openxlsx)
library(broom)
library(psych)
library(emmeans)
library(knitr)
library(rstatix)
library(performance)
library(ggseg)
library(grid)
```

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

```

# 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, {
  education_c <- C(as.factor(education), treatment, 3)
  print(attributes(education_c))
})
  
```

## Thickness and cognition
```{r, fig.height=9, fig.width=8}
## Thickness data
all_data_clean[140] <- read.table("data/frag_index_sd.csv", col.names = "fi_thickness_sd")
all_data_clean[141] <- read.table("data/frag_index_sd_avg.csv", col.names = "fi_thickness_sd_avg")

## Select variable names for models
outcome_var <- c("adas_cog_tot","dsst", "tmt_bma_lg")

# List to store models
cortex_lm = list()

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

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

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

## Sleep and cognition
### Sleep duration
#### Average

```{r, fig.height=9, fig.width=8}
# List to store models
sleep_duration_avg_lm = list()

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

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

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

```

#### Variability
```{r, fig.height=9, fig.width=8}
# List to store models
sleep_duration_sd_lm = list()

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

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

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

### Sleep efficiency
#### Average

```{r, fig.height=9, fig.width=8}
# List to store models
sleep_efficiency_avg_lm = list()

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

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

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

#### Variability

```{r, fig.height=9, fig.width=8}
# List to store models
sleep_efficiency_sd_lm = list()

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

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

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

### Sleep fragmentation
#### Average

```{r, fig.height=9, fig.width=8}
# List to store models
fragmentation_index_avg_lm = list()

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

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

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

#### Variability

```{r, fig.height=9, fig.width=8}
# List to store models
fragmentation_index_sd_lm = list()

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

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

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

## Tables
### Thickness models

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

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

# Print estimates
cortex_lm_estimates %>% 
  kable(digits = 2)
```

### Sleep models

```{r}
## List with all models
sleep_lm <- c(sleep_duration_avg_lm, sleep_duration_sd_lm, 
              sleep_efficiency_avg_lm, sleep_efficiency_sd_lm, 
              fragmentation_index_avg_lm, fragmentation_index_sd_lm)
## Adjusted
sleep_lm_estimates = tibble()

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

# Print estimates
sleep_lm_estimates %>% 
  kable(digits = 2)
```

## Significant clusters

```{r, fig.height=9, fig.width=8}
### Superior frontal gyrus
### Model C: Age, sex, MoCA, Sleep Medication

### Fragmentation index variability
lm_fragmnt_index_model_c <- lm(fi_thickness_sd ~ fragmentation_index_sd_lg
                                  + age
                                  + sex
                                  + moca
                                  + sleep_medication_psqi, all_data_clean)
lm_results(lm_fragmnt_index_model_c)
check_model(lm_fragmnt_index_model_c)

### Fragmentation index variability
lm_fragmnt_index_model_d <- lm(fi_thickness_sd_avg ~ fragmentation_index_sd_lg
                                  + age
                                  + sex
                                  + moca
                                  + sleep_medication_psqi
                                  + fragmentation_index_avg, all_data_clean)
lm_results(lm_fragmnt_index_model_d)
check_model(lm_fragmnt_index_model_d)
```

### Figures

```{r }
### Model c
ggplot(all_data_clean,aes(fragmentation_index_sd_lg, fi_thickness_sd)) +
  geom_point(size = 3, colour = "black") +
  ggtitle("Left superior frontal gyrus") +
  geom_smooth(method = lm, se = TRUE, colour = "black", fill = "green") +
  xlab("Fragmentation index variability (log transformed)") +
  ylab(bquote('Thickness,'~mm^2)) + 
  scale_y_continuous(breaks = seq(2, 3.5, by = 0.25), limits = c(2, 3.25)) +
  theme_light() +
  theme(plot.title = element_text(size = 17, hjust = 0.5),
        axis.line = element_blank(),
        axis.title.y = element_text (size = 15),
        axis.title.x = element_text (size = 15),
        axis.text = element_text(size = 13))
ggsave("plots/model_a.png", width = 7, height = 5, dpi = 300)

### Model d
ggplot(all_data_clean,aes(fragmentation_index_sd_lg, fi_thickness_sd_avg)) +
  geom_point(size = 3, colour = "black") +
  ggtitle("Left superior frontal gyrus") +
  geom_smooth(method = lm, se = TRUE, colour = "black", fill = "turquoise") +
  xlab("Fragmentation index variability (log transformed)") +
  ylab(bquote('Thickness,'~mm^2)) + 
  scale_y_continuous(breaks = seq(2, 3.5, by = 0.25), limits = c(2, 3.25)) +
  theme_light() +
  theme(plot.title = element_text(size = 17, hjust = 0.5),
        axis.line = element_blank(),
        axis.title.y = element_text (size = 15),
        axis.title.x = element_text (size = 15),
        axis.text = element_text(size = 13))
ggsave("plots/model_b.png", width = 7, height = 5, dpi = 300)
```


# Moderation analysis
```{r}
# Computing mean thickness for each roi
thickness_data <- all_data_clean %>% 
  select(id, age,sex, moca, sleep_medication_psqi, education_c, adas_cog_tot, dsst, tmt_bma_lg, 
         fragmentation_index_avg, fragmentation_index_sd_lg, 
         contains("thickness"), 
         -contains("mean_thickness"), -fi_thickness_sd, -fi_thickness_sd_avg, -contains("model")) %>%
  pivot_longer(contains("thickness")) %>% 
  mutate(hemisphere = ifelse(str_detect(name, "lh_"), "left", "right"),
         name = str_remove(name, "lh_|rh_")) %>% 
  pivot_wider(names_from = hemisphere, values_from = value) %>% 
  mutate(mean_thickness = (left + right) / 2) %>% 
  select(-left, -right) %>% 
  pivot_wider(values_from = mean_thickness)

# ROIs
rois <- names(thickness_data %>% select(contains("thickness")))
```

## ADAS-Cog
### Average

```{r}
# List to store models
thickness_adas_avg_lm = list()

thickness_data <- thickness_data %>% 
  mutate(education_c2 = ifelse(education_c == "university degree", "university degree", "no university"))

# Run models in loop script
for (var in rois){
  
  # Formula
  form <- paste(var,"~ fragmentation_index_sd_lg + age + sex + education_c2 + sleep_medication_psqi")

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

## Variability
for (model in thickness_adas_avg_lm){
  print(formula(model), showEnv = FALSE)
  print(lm_results(model))
}

summary(thickness_data)
```

### Variability

```{r}
# List to store models
thickness_adas_sd_lm = list()

# Run models in loop script
for (var in rois){
  
  # Formula
  form <- paste("adas_cog_tot","~",var,"* fragmentation_index_sd_lg + age + sex + sleep_medication_psqi + education_c")

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

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

## DSST
### Average

```{r}
# List to store models
thickness_dsst_avg_lm = list()

# Run models in loop script
for (var in rois){
  
  # Formula
  form <- paste("dsst","~",var,"* fragmentation_index_avg + age + sex + sleep_medication_psqi + education_c")

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

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

### Variability

```{r}
# List to store models
thickness_dsst_sd_lm = list()

# Run models in loop script
for (var in rois){
  
  # Formula
  form <- paste("dsst","~",var,"* fragmentation_index_sd_lg + age + sex + sleep_medication_psqi + education_c")

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

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

## TMT
### Average

```{r}
# List to store models
thickness_tmt_avg_lm = list()

# Run models in loop script
for (var in rois){
  
  # Formula
  form <- paste("tmt_bma_lg","~",var,"* fragmentation_index_avg + age + sex + sleep_medication_psqi + education_c")

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

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

### Variability

```{r}
# List to store models
thickness_tmt_sd_lm = list()

# Run models in loop script
for (var in rois){
  
  # Formula
  form <- paste("tmt_bma_lg","~",var,"* fragmentation_index_sd_lg + age + sex + sleep_medication_psqi + education_c")

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

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

\newpage
\blandscape

## Tables

```{r} 
# Average
cortex_avg_mod = list()
cortex_avg_mod <- c(thickness_adas_avg_lm,
                    thickness_dsst_avg_lm,
                    thickness_tmt_avg_lm)

cortex_avg_mod_estimates = tibble()

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

## Compute FDR by outcome (tests = 34)
cortex_avg_mod_estimates <- cortex_avg_mod_estimates %>% 
  separate(term, c("predictor", "moderator"), ":") %>% 
  group_by(outcome) %>% 
  mutate(p.fdr = p.adjust(p.value, method = "fdr")) %>% 
  ungroup()

# Print models with FDR p values
cortex_avg_mod_estimates %>% 
  filter(p.value <= 0.05) %>% 
  arrange(outcome, p.value) %>% 
  kable(digits = 3)

# Variability
cortex_sd_mod = list()
cortex_sd_mod <- c(thickness_adas_sd_lm,
                  thickness_dsst_sd_lm,
                  thickness_tmt_sd_lm)

cortex_sd_mod_estimates = tibble()

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

## Compute FDR by outcome (tests = 34)
cortex_sd_mod_estimates <- cortex_sd_mod_estimates %>% 
  separate(term, c("predictor", "moderator"), ":") %>% 
  group_by(outcome) %>% 
  mutate(p.fdr = p.adjust(p.value, method = "fdr")) %>% 
  ungroup()

# Print models with FDR p values
cortex_sd_mod_estimates %>% 
  filter(p.value <= 0.05) %>% 
  arrange(outcome, p.value) %>% 
  kable(digits = 2)

rbind(cortex_avg_mod_estimates, cortex_sd_mod_estimates) %>% 
  filter(p.value <= 0.05) %>% 
  arrange(predictor) %>%
  distinct(predictor) %>% 
  mutate(index = rank(predictor)) %>% 
  kable()
```

\elandscape
\newpage

## Figures
### Cortex

```{r, fig.height=10, fig.width=6}
# Average
## Add labels for lh and rh 
plot_cortex_avg <- rbind(
  cortex_avg_mod_estimates %>% 
  mutate(predictor = paste("lh_", predictor, sep = ""),
         predictor = str_remove(predictor, "_thickness")),
  
  cortex_avg_mod_estimates %>% 
  mutate(predictor = paste("rh_", predictor, sep = ""),
         predictor = str_remove(predictor, "_thickness"))
)

## Plot significant regions
plot_avg <- plot_cortex_avg %>% 
  mutate(t.critical = ifelse(statistic >= 1.96, statistic, ifelse(statistic <= -1.96, statistic, 0))) %>% 
  group_by(outcome) %>% 
  rename(label = predictor) %>% 
  mutate(outcome = str_replace_all(
    outcome, c("adas_cog_tot" = "ADAS-Cog",
               "dsst" = "DSST",
               "tmt_bma_lg" = "TMT (B - A)"))) %>% 
  ggplot() +
  geom_brain(atlas = dk, 
             aes(fill = t.critical),
              position = position_brain(side ~ hemi), show.legend = TRUE) + 
  scale_fill_gradient2(
    high = "#FFCF00", low = "#8C6DD5", mid = "grey90", name = "Uncorrected p\n",
    midpoint = 0,
    labels = c("0", "0.01", "0.05", "0.05", "0.01"),
    breaks = c(0, -2.576, -1.96, 1.96, 2.576) , limits = c(-3, 3)) +
  facet_wrap(~outcome, ncol = 1) +
  ggtitle("Fragmentation Index Average") +
  theme_minimal() +
  theme(plot.background = element_rect(fill = "white", colour = NA),
        plot.title = element_text(hjust = 0.5, size = 10),
        axis.text = element_blank(),
        legend.position = "none")

# Variability
plot_cortex_sd <- rbind(
  cortex_sd_mod_estimates %>% 
  mutate(predictor = paste("lh_", predictor, sep = ""),
         predictor = str_remove(predictor, "_thickness")),
  
  cortex_sd_mod_estimates %>% 
  mutate(predictor = paste("rh_", predictor, sep = ""),
         predictor = str_remove(predictor, "_thickness"))
)

## Add labels for lh and rh 
plot_sd <-plot_cortex_sd %>% 
  mutate(t.critical = ifelse(statistic >= 1.96, statistic, ifelse(statistic <= -1.96, statistic, 0))) %>% 
  group_by(outcome) %>% 
  rename(label = predictor) %>% 
  mutate(outcome = str_replace_all(
    outcome, c("adas_cog_tot" = "ADAS-Cog",
               "dsst" = "DSST",
               "tmt_bma_lg" = "TMT (B - A)"))) %>% 
  ggplot() +
  geom_brain(atlas = dk, 
             aes(fill = t.critical),
              position = position_brain(side ~ hemi), show.legend = TRUE) +
  scale_fill_gradient2(
    high = "#FFCF00", low = "#8C6DD5", mid = "grey90", name = "Uncorrected p\n",
    midpoint = 0,
    labels = c("0", "Negative", "0.01", "0.05", "0.05", "0.01", "Postive"),
    breaks = c(0, -3, -2.576, -1.96, 1.96, 2.576, 3) , limits = c(-3, 3)) +
  facet_wrap(~outcome, ncol = 1) +
  ggtitle("Fragmentation Index Variability") +
  theme_minimal() +
  theme(plot.background = element_rect(fill = "white", colour = NA),
        plot.title = element_text(hjust = 0.5, size = 10),
        axis.text = element_blank(), 
        legend.title = element_text(size = 7),
        legend.text = element_text(size = 6),
        legend.position = "none",
        legend.key.width = unit(0.3, "cm"))

plot_avg + plot_sd

# Merging and saving plots
ggsave(plot = (plot_avg + plot_sd), width = 7, height = 8, "plots/thickness_mod.png")
```

### Scatterplots

```{r}
# Creating grouping variables to plot moderation based on above or below median
thickness_data <- thickness_data %>% 
  mutate(group_sd = ifelse(fragmentation_index_sd_lg > median(fragmentation_index_sd_lg), "High", "Low"),
         group_avg = ifelse(fragmentation_index_avg > median(fragmentation_index_avg), "High", "Low"))
```

```{r, fig.height=4, fig.width=4}
# Plotting average models
thickness_data %>% 
  ggplot(aes(precentral_thickness, dsst)) +
  geom_point(size = 1.5) +
  geom_smooth(aes(colour = group_avg), method = "lm") +
  ggtitle("Precentral gyrus") +
  labs(y = "DSST", x = "Cortical thickness") +
  scale_colour_viridis_d(name = "Fragmentaition index \naverage") +
  theme_minimal() +
  theme(plot.background = element_rect(fill = "white", colour = NA),
        plot.title = element_text(hjust = 0.5, size = 10),
      legend.position = "bottom",
      axis.title = element_text (size = 10),
      axis.text = element_text(size = 6))
ggsave("plots/dsst_avg.png", dpi = 300, height = 3, width = 4)
```

```{r, fig.height=6, fig.width=4}
thickness_data %>%
  pivot_longer(c(posteriorcingulate_thickness, superiorfrontal_thickness)) %>% 
  mutate(name = str_remove(str_to_sentence(name), "_thickness")) %>% 
  ggplot(aes(value, tmt_bma_lg)) +
  geom_point(size = 1.5) +
  geom_smooth(aes(colour = group_avg), fill = "gray90", method = "lm") +
  scale_colour_viridis_d(name = "Fragmentaition index \naverage") +
  labs(y = "TMT (B - A)", x = "Cortical thickness") +
  facet_wrap(~name, ncol = 1, scales = "free_x") +
  theme_minimal() +
  theme(plot.background = element_rect(fill = "white", colour = NA),
        plot.title = element_text(hjust = 0.5, size = 10),
      legend.title = element_text(size = 6),
      legend.text = element_text(size = 6),
      legend.position = "bottom",
      axis.title = element_text (size = 10),
      axis.text = element_text(size = 6))
ggsave("plots/tmt_avg.png", dpi = 300, height = 6, width = 4)

```

```{r, fig.height=4, fig.width=4}
# Variability
# Plotting variability models
thickness_data %>% 
  ggplot(aes(bankssts_thickness, adas_cog_tot)) +
  geom_point(size = 1.5) +
  geom_smooth(aes(colour = group_sd), method = "lm") +
  scale_colour_viridis_d() +
  ggtitle("Banks of the superior temporal sulcus") +
  labs(y = "ADAS-Cog", x = "Cortical thickness") +
  theme_minimal() +
  theme(plot.background = element_rect(fill = "white", colour = NA),
        plot.title = element_text(hjust = 0.5, size = 10),
      legend.title = element_text(size = 6),
      legend.text = element_text(size = 6),
      legend.position = "bottom",
      axis.title = element_text (size = 10),
      axis.text = element_text(size = 6))
ggsave("plots/adas_sd.png", dpi = 300, height = 3, width = 4)
```

```{r, fig.height=10, fig.width=4}
thickness_data %>%
  pivot_longer(c(cuneus_thickness, pericalcarine_thickness, precentral_thickness, precuneus_thickness)) %>% 
  mutate(name = str_remove(str_to_sentence(name), "_thickness")) %>% 
  ggplot(aes(value, dsst)) +
  geom_point(size = 1.5) +
  labs(y = "DSST", x = "Cortical thickness") +
  geom_smooth(aes(colour = group_sd), fill = "gray90", method = "lm") +
  scale_colour_viridis_d(name = "Fragmentaition index \navariability") +
  facet_wrap(~name, ncol = 1, scales = "free_x") +
  theme_minimal() +
  theme(plot.background = element_rect(fill = "white", colour = NA),
        plot.title = element_text(hjust = 0.5, size = 10),
      legend.title = element_text(size = 6),
      legend.text = element_text(size = 6),
      legend.position = "bottom",
      axis.title = element_text (size = 10),
      axis.text = element_text(size = 6))
ggsave("plots/dsst_sd.png", dpi = 300, height = 12, width = 4)
```

```{r, fig.height=10, fig.width=4}

thickness_data %>%
  pivot_longer(c(cuneus_thickness, fusiform_thickness, supramarginal_thickness, transversetemporal_thickness)) %>% 
  mutate(name = str_remove(str_to_sentence(name), "_thickness")) %>% 
  ggplot(aes(value, dsst)) +
  geom_point(size = 1.5) +
  labs(y = "TMT (B - A)", x = "Cortical thickness") + 
  geom_smooth(aes(colour = group_sd), fill = "gray90", method = "lm") +
  scale_colour_viridis_d(name = "Fragmentaition index avariability") +
  facet_wrap(~name, ncol = 1, scales = "free_x") +
  theme_minimal() +
 theme(plot.background = element_rect(fill = "white", colour = NA),
      plot.title = element_text(hjust = 0.5, size = 10),
      legend.title = element_text(size = 6),
      legend.text = element_text(size = 6),
      legend.position = "bottom",
      axis.title = element_text (size = 10),
      axis.text = element_text(size = 6))
ggsave("plots/tmt_sd.png", dpi = 300, height = 12, width = 4)
```