ComBat_application.Rmd

---
title: "Supplementary information - COMBAT testing"
author: "Fernanda Hansen Pacheco de Moraes"
date: "18 jan 2022"
output:
  html_document: 
    fig_caption: yes
    fig_width: 8
    number_sections: yes
    theme: paper
    toc: yes
editor_options: 
  chunk_output_type: inline
---

Description of the procedures and analysis present in Manuscript 2,
**Establishing a baseline for human cortical folding morphological variables: a multicenter study**, at the Doctorate Thesis presented
to the Programa de Pós-Graduação em Ciências Médicas at the Instituto
D'Or de Pesquisa e Ensino as a partial requirement to obtain the
Doctorate Degree.

Part of the data used here cannot be shared due to restrictions of the
Ethic Committee. Data can be shared upon reasonable request to the
corresponding author. To fulfill these limitation, we will generate
random data to simulate the results.

Get in touch with us
([fernandahmoraes\@gmail.com](mailto:fernandahmoraes@gmail.com){.email})
in case any help is needed, our aim is to improve the code as needed!

# RMARKDOWN AND R SETUP
```{r setup, include=FALSE}
knitr::opts_chunk$set(
	echo = TRUE,
	message = FALSE,
	warning = FALSE,
	cache = TRUE
)
```

```{r working directory}
setwd("D:/GitHub/Typical-values")
```

```{r functions, message=FALSE, warning=FALSE}
## define functions

# test angular coeficinet versus theoretical value
test_coef <- function(reg, coefnum, val){
  co <- coef(summary(reg))
  tstat <- (co[coefnum,1] - val)/co[coefnum,2]
  2 * pt(abs(tstat), reg$df.residual, lower.tail = FALSE)
}

# wrap text
wrapper <- function(x, ...) paste(strwrap(x, ...), collapse = "\n")
```

```{r call packages}
library(readr)
library(tidyverse)
library(lubridate)
library(ggpubr)
library(kableExtra)
library(broom)
library(lme4)
library(lsmeans)
library(MuMIn)
library(arm)
library(effects)
library(compute.es)

library(devtools)
# install_github("jfortin1/neuroCombatData")
# install_github("jfortin1/neuroCombat_Rpackage")
library(neuroCombat)

```

```{r} 
# COLOR BLIND PALETTE WITH BLACK
cbbPalette <- c("#000000", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
cbbPalette2 <- c("#56B4E9", "#E69F00", "#CC79A7", "#56B4E9", "#D55E00", "#0072B2", "#009E73", "#F0E442")
```

# set seed for random process

```{r}
set.seed(1)
```

```{r import files}
# dados_datasetscomp <- read_csv("dados_datasetscomp2.csv")
dados_datasetscomp <- read_csv("data_typical_values.csv")

dados_datasetscomp <- dados_datasetscomp %>%
  filter(Sample != "HCP500r", Diagnostic == "CTL", ROI == "hemisphere")

```

# DATA PREPARATION
```{r create new variables}
# estimate cortical folding variables
dados_datasetscomp <- dados_datasetscomp %>%
  mutate(
    # create new variables
    GMvolume = ifelse(!is.na(GMvolume),GMvolume,AvgThickness*TotalArea),
    logAvgThickness = log10(AvgThickness),
    logTotalArea = log10(TotalArea),
    logExposedArea = log10(ExposedArea),
    localGI = TotalArea / ExposedArea,
    k = sqrt(AvgThickness) * TotalArea / (ExposedArea ^ 1.25),
    K = 1 / 4 * log10(AvgThickness ^ 2)  + log10(TotalArea) - 5 / 4 * log10(ExposedArea),
    S = 3 / 2 * log10(TotalArea) + 3 / 4 * log10(ExposedArea) - 9 /  4 * log10(AvgThickness ^
                                                                                 2) ,
    I = log10(TotalArea) + log10(ExposedArea) + log10(AvgThickness ^ 2),
    # c = as.double(ifelse(
    #   ROI == "hemisphere", NA, 4 * pi / GaussianCurvature
    # )),
    TotalArea_corrected = ifelse(ROI == "hemisphere", TotalArea, TotalArea * c),
    ExposedArea_corrected = ifelse(ROI == "hemisphere", ExposedArea, ExposedArea * c),
    logTotalArea_corrected = log10(TotalArea_corrected),
    logExposedArea_corrected = log10(ExposedArea_corrected),
    localGI_corrected = ifelse(
      ROI == "hemisphere",
      TotalArea / ExposedArea,
      TotalArea_corrected / ExposedArea_corrected
    ),
    k_corrected = ifelse(
      ROI == "hemisphere",
      sqrt(AvgThickness) * log10(TotalArea) / (log10(ExposedArea) ^ 1.25),
      sqrt(AvgThickness) * log10(TotalArea_corrected) / (log10(ExposedArea_corrected ^
                                                                 1.25))
    ),
    K_corrected =  ifelse(
      ROI == "hemisphere",
      1 / 4 * log10(AvgThickness ^ 2) + log10(TotalArea) - 5 / 4 * log10(ExposedArea),
      1 / 4 * log10(AvgThickness ^ 2) + log10(TotalArea_corrected) - 5 / 4 * log10(ExposedArea_corrected)
    ),
    I_corrected = ifelse(
      ROI == "hemisphere",
      log10(TotalArea) + log10(ExposedArea) + log10(AvgThickness ^ 2) ,
      log10(TotalArea_corrected) + log10(ExposedArea_corrected) + log10(AvgThickness ^ 2)
    ),
    S_corrected = ifelse(
      ROI == "hemisphere",
      3 / 2 * log10(TotalArea) + 3 / 4 * log10(ExposedArea) - 9 /  4 * log10(AvgThickness ^ 2) ,
      3 / 2 * log10(TotalArea_corrected) + 3 / 4 * log10(ExposedArea_corrected) - 9 /  4 * log10(AvgThickness ^ 2)
    ),
    Knorm = K_corrected / sqrt(1 + (1 / 4) ^ 2 + (5 / 4) ^ 2),
    Snorm = S_corrected / sqrt((3 / 2) ^ 2 + (3 / 4) ^ 2 + (9 / 4) ^ 2),
    Inorm = I_corrected / sqrt(1 ^ 2 + 1 ^ 2 + 1 ^ 1)
  )

# create age intervals
dados_datasetscomp$Age_interval <- cut(dados_datasetscomp$Age,
                                       breaks = c(0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, 100),
                                       right = FALSE,
                                       include.lowest = TRUE)

dados_datasetscomp$Age_interval10 <- cut(dados_datasetscomp$Age,
                                         breaks = c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100),
                                         right = FALSE,
                                         include.lowest = TRUE)
```

```{r data cleanup}
dados_all <- dados_datasetscomp %>% filter(!is.na(logAvgThickness), ExposedArea != 0 | !is.na(localGI), !is.infinite(logExposedArea)) %>% 
  droplevels()

dados_datasetscomp <- dados_all

dados_datasetscomp$Diagnostic <- as.factor(dados_datasetscomp$Diagnostic)
dados_datasetscomp$Diagnostic <- relevel(dados_datasetscomp$Diagnostic, ref = "CTL")

```

# Comparing with our harmonization

## AvgThickness

Here we included hemisphere as a factor since the combat needs two lines of morphological measures input.
```{r}

m.1 <-
  lme4::lmer(AvgThickness ~ Age + hemi + (1|Sample), data = dados_datasetscomp)

equatiomatic::extract_eq(m.1, wrap=TRUE)

fixef(m.1)

re <- as_tibble(ranef(m.1)) %>%
  filter(grpvar == "Sample") %>%
  mutate(
    T_shift = condval,
    sd_T_shift = condsd,
    Sample = grp
  ) %>%
  dplyr::select(-c(condval, grpvar, term, condsd, grp))

dados_datasetscomp_harmonized <- dados_datasetscomp %>%
   full_join(re) %>%
   mutate(
     AvgThickness_shiftc = AvgThickness - T_shift
   )

```

```{r}
summary(m.1)
sigma(m.1)
sigma.hat(m.1)
```

```{r}
var(dados_datasetscomp$AvgThickness, dados_datasetscomp$Age)
```


```{r}
re


b <- ggplot(re, aes(Sample, T_shift)) +
  geom_point() +
  geom_errorbar(aes(ymin = T_shift - sd_T_shift, ymax  = T_shift + sd_T_shift)) +
  theme_pubr() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

b
```

## After hamonization

```{r}
m.1 <-
  lme4::lmer(AvgThickness_shiftc ~ Age + hemi + (1|Sample) , data = dados_datasetscomp_harmonized)

equatiomatic::extract_eq(m.1, wrap=TRUE)

fixef(m.1)

re <- as_tibble(ranef(m.1)) %>%
  filter(grpvar == "Sample") %>%
  mutate(
    T_shift = condval,
    sd_T_shift = condsd,
    Sample = grp
  ) %>%
  dplyr::select(-c(condval, grpvar, term, condsd, grp))

summary(m.1)
sigma(m.1)
sigma.hat(m.1)

var(dados_datasetscomp_harmonized$AvgThickness_shiftc, dados_datasetscomp_harmonized$Age)

as_tibble(ranef(m.1)) %>%
  filter(grpvar == "hemi")

d <- ggplot(re, aes(Sample, T_shift)) +
  geom_point() +
  geom_errorbar(aes(ymin = T_shift - sd_T_shift, ymax  = T_shift + sd_T_shift)) +
  theme_pubr() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  ylim(-0.05, 0.1)

d
```

```{r}
a <- ggplot(dados_datasetscomp_harmonized, aes(x = Age, color = Sample, alpha = 0.4)) + 
  geom_point(aes(y = AvgThickness)) + 
  guides(alpha = "none") +
  theme_pubr()

a

c <- ggplot(dados_datasetscomp_harmonized, aes(x = Age, y = AvgThickness_shiftc, color = Sample, alpha = 0.4)) + 
  geom_point() + 
  guides(alpha = "none") +
  theme_pubr()

c
```

# Combat preparation

```{r}
## Here we will select only cortical thickness for hemispheres

dados_datasetscomp <- dados_datasetscomp %>%
  filter(ROI == "hemisphere") %>%
  dplyr::select(c(SUBJ, hemi, Sample, Age, AvgThickness)) %>%
  pivot_wider(names_from = hemi, values_from = AvgThickness)

dat <- as.data.frame(t(as.matrix(dplyr::select(dados_datasetscomp, c(L, R)))))

batch <- as.data.frame(t(as.matrix(dplyr::select(dados_datasetscomp, c(Sample)))))

mod <- model.matrix(~Age, dados_datasetscomp)

```

```{r}
results <- neuroCombat(dat, batch, mod)
```

# COMBAT RESULTS

```{r}
as.data.frame(results$estimates$gamma.hat)
as.data.frame(results$estimates$var.pooled)
```

```{r}
neuroCombat::drawPriorDelta(results$estimates)
neuroCombat::drawPriorGamma(results$estimates)
```

```{r}
data_standardized <- results[["dat.standardized"]]
col <- as.factor(batch)
boxplot(data_standardized, ylim=c(-2,2),col=as.factor(batch))
```


```{r}
data_norm <- results[["dat.combat"]]
data_norm <- standardizeData(data_norm,
                             batch=batch,
                             mod=mod)
boxplot(data_norm, ylim=c(-2,2),col=col)
```

## After hamonization

```{r}
AvgThickness_shiftc_combat <- as_tibble(t(results$dat.combat)) %>% 
  cbind(mod) %>%
  cbind(as_tibble(t(results$estimates$batch))) %>%
  pivot_longer(c(L, R), names_to = "hemi", values_to = "AvgThickness_shiftc_combat") %>%
  dplyr::select(-c(`(Intercept)`))

AvgThickness_combat  <- as_tibble(t(results$dat.original)) %>% 
  cbind(mod) %>%
  cbind(as_tibble(t(results$estimates$batch))) %>%
  pivot_longer(c(L, R), names_to = "hemi", values_to = "AvgThickness_combat") %>%
  dplyr::select(-c(`(Intercept)`))

harmonized_combat <- full_join(AvgThickness_combat, AvgThickness_shiftc_combat)

```
### confirming the results before harmonization, with the data from combat product

```{r}
m.1 <-
  lme4::lmer(AvgThickness_combat ~ Age + hemi + (1|Sample) , data = harmonized_combat)

equatiomatic::extract_eq(m.1, wrap=TRUE)

fixef(m.1)

re <- as_tibble(ranef(m.1)) %>%
  filter(grpvar == "Sample") %>%
  mutate(
    T_shift = condval,
    sd_T_shift = condsd,
    Sample = grp
  ) %>%
  dplyr::select(-c(condval, grpvar, term, condsd, grp))

summary(m.1)
sigma(m.1)
sigma.hat(m.1)

var(harmonized_combat$AvgThickness_combat, harmonized_combat$Age)

as_tibble(ranef(m.1)) %>%
  filter(grpvar == "hemi")


b_c <- ggplot(re, aes(Sample, T_shift)) +
  geom_point() +
  geom_errorbar(aes(ymin = T_shift - sd_T_shift, ymax  = T_shift + sd_T_shift)) +
  theme_pubr() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

b_c
```

#### and then the harmonized data

```{r}
m.1 <-
  lme4::lmer(AvgThickness_shiftc_combat ~ Age + hemi + (1|Sample) , data = harmonized_combat)

equatiomatic::extract_eq(m.1, wrap=TRUE)

fixef(m.1)

re <- as_tibble(ranef(m.1)) %>%
  filter(grpvar == "Sample") %>%
  mutate(
    T_shift = condval,
    sd_T_shift = condsd,
    Sample = grp
  ) %>%
  dplyr::select(-c(condval, grpvar, term, condsd, grp))

summary(m.1)
sigma(m.1)
sigma.hat(m.1)

var(harmonized_combat$AvgThickness_shiftc_combat, harmonized_combat$Age)

as_tibble(ranef(m.1)) %>%
  filter(grpvar == "hemi")

d_c <- ggplot(re, aes(Sample, T_shift)) +
  geom_point() +
  geom_errorbar(aes(ymin = T_shift - sd_T_shift, ymax  = T_shift + sd_T_shift)) +
  theme_pubr() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  ylim(-0.05, 0.1)

d_c
```

```{r}
a_c <- ggplot(harmonized_combat, aes(x = Age, color = Sample, alpha = 0.4)) + 
  geom_point(aes(y = AvgThickness_combat)) +
  guides(alpha = "none") +
  theme_pubr()

a_c

c_c <- ggplot(harmonized_combat, aes(x = Age, color = Sample, alpha = 0.4)) + 
  geom_point(aes(y = AvgThickness_shiftc_combat)) + 
  guides(alpha = "none") +
  theme_pubr()

c_c
```

### fazer figura com os ranef e com antes e depois tudo junto
```{r}
fig_combat_1_1 <- ggarrange(b, d, labels = c("A", "B"), nrow = 2, ncol = 1, font.label = list(size = 11), common.legend = TRUE, legend = "top", align = "hv") %>%
  annotate_figure(top = text_grob("Our hamonization"))

fig_combat_1_2 <- ggarrange(b_c, d_c, labels = c("C", "D"), nrow = 2, ncol = 1, font.label = list(size = 11), common.legend = TRUE, legend = "top", align = "hv") %>%
  annotate_figure(top = text_grob("Combat"))

fig_combat_1 <- ggarrange(fig_combat_1_1, fig_combat_1_2, ncol = 2, font.label = list(size = 11), common.legend = TRUE, legend = "top", align = "hv") 

fig_combat_1

ggsave("fig_combat_1.pdf", plot = fig_combat_1, width = 18, height = 22, units = "cm", device = "pdf")
```

```{r}
fig_combat_2_1 <- ggarrange(a, c, labels = c("A", "B"), nrow = 1, ncol = 2, font.label = list(size = 11), common.legend = TRUE, legend = "top", align = "hv") %>%
  annotate_figure(top = text_grob("Our hamonization"))

fig_combat_2_2 <- ggarrange(a_c, c_c, labels = c("C", "D"), nrow = 1, ncol = 2, font.label = list(size = 11), common.legend = TRUE, legend = "top", align = "hv") %>%
  annotate_figure(top = text_grob("Combat"))

fig_combat_2 <- ggarrange(fig_combat_2_1, fig_combat_2_2, ncol = 1, font.label = list(size = 11), common.legend = TRUE, legend = "top", align = "hv") 

fig_combat_2

ggsave("fig_combat_2.png", plot = fig_combat_2, width = 18, height = 22, units = "cm", device = "png")
ggsave("fig_combat_2.pdf", plot = fig_combat_2, width = 18, height = 22, units = "cm", device = "pdf")
```