---
title: "BrainAge_VidalPineiro_Tests(Lifebrain replication)"
author: "dvp"
date: "11/2/2020"
output:
  html_document: default
  pdf_document: default
---

# Setup
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(tidyverse)
library(magrittr)
#library(itsadug)
#library(lavaan)
#library(reshape2)
library(broom)
library(lmerTest)
library(MuMIn)


data_folder=file.path("./data/LifeBrain",paste0("noExcl_resCohort_outlierremove"))
options(bitmapType = "cairo")
```

# Loading data and common Preproc
load data
```{r load}
load(file.path(data_folder, "ResultsBrainAge.Rda"))
df.Test = df.Test[[1]]
df.Test %<>% mutate(BAc = age + BAG_corr_gam)
```

# Compute change
Delta considered as the GAM-coefficient corrected GAP
Change in Delta derived from linear models

```{r compute delta change}

# compute interceptect (demeaned) and slope
tmp.BAG_corr = df.Test %>% 
  group_by(CrossProject_ID) %>% 
  mutate(Time = Time - mean(Time)) %>% 
  do(fit = lm(BAG_corr_gam ~ Time, data = .))

tmp.BAG_corr$df = lapply(tmp.BAG_corr$fit, tidy)  

df.lm = data.frame(CrossProject_ID = tmp.BAG_corr$CrossProject_ID,
                BAG_corr_mean = lapply(tmp.BAG_corr$df, 
                                      function(x) x$estimate[[1]]) %>% simplify2array(),
                BAG_corr_change = lapply(tmp.BAG_corr$df, 
                                         function(x) x$estimate[[2]]) %>% simplify2array())


df.lm2 = df.Test %>% 
  group_by(CrossProject_ID) %>% 
  summarise(sex = first(sex),
             AgeBsl = first(AgeBsl),
             cohort = first(cohort),
            eICV = mean(Vol_EstimatedTotalIntraCranial_wb),
            AgeGap = first(AgeGap),
            n = n())

df.Test.long = left_join(df.lm, df.lm2)
```

# Analysis
### a) Cross/sectional delta predicts brain aging
```{r Brain aging to aging}
## BAG (Brain Age Gap bias-corrected at tp1 predicts less Brain Aging Change with time)
lm.pred = lm(BAG_corr_change ~ BAG_corr_mean, data = df.Test.long) # uncorrected. just to check 
summary(lm.pred)

# as in UKB
lm.main = lmer(BAG_corr_change ~ BAG_corr_mean + AgeBsl + sex + scale(eICV) + (1 | cohort), data = df.Test.long, REML = T)
summary(lm.main)

  # save  residuals
  df.Test.long$residuals = summary(lm.main)$residuals %>% as.numeric()
  
  # UVE - unique variance explained
  lm.main_rm = lmer(BAG_corr_change ~  AgeBsl + sex + scale(eICV) + (1 | cohort), data = df.Test.long, REML = T)
  r.squaredGLMM(lm.main) -r.squaredGLMM(lm.main_rm)
  
  # visually check assumptions
  plot(lm.main)
  qqnorm(residuals(lm.main))
  plot(lm.main@frame$BAG_corr_mean,residuals(lm.main))

# supplementary analysis controling for age gap
# control for agegap and restrict analysis to 4 years gap
lm.main_agegap1 = lmer(BAG_corr_change ~ BAG_corr_mean + AgeBsl + AgeGap + sex + scale(eICV) + (1 | cohort), data = df.Test.long, REML = T)
summary(lm.main_agegap1)

sum(df.lm2$AgeGap > 4)
lm.main_agegap2 = lmer(BAG_corr_change ~ BAG_corr_mean + AgeBsl + sex + scale(eICV) + (1 | cohort), data = df.Test.long %>% filter(AgeGap > 4), REML = T)
summary(lm.main_agegap2)
```

### b) Fig. 2b. Plotting cross to long
```{r Plot cross to long}
# Lifebrain - XGB
gs = ggplot(data = df.Test.long,
       mapping = aes(x = BAG_corr_mean,
                     y = residuals)) + 
  geom_point(shape = 21, color = "black", fill = "#66a61e", size = 4, stroke = 1.5) + 
  geom_smooth(method = "lm", 
              color = "#1b9e77",
              size = 3) +
  theme_classic() + 
  theme(legend.position = 'none',
        axis.text = element_text(size = 16),
        axis.title = element_text(size = 20),
        plot.title = element_text(size = 20, hjust = 0.5, vjust = 0)) +
  ylab(expression("Brain Age Delta"[long])) +
  xlab(expression("Brain Age Delta"[cross])) +
  ggtitle("Lifebrain - XGB")

ggsave("figures/BrainAge_Lifebrain.png", 
       dpi = 500, 
       plot = gs,
       width = 10,
       height = 10, 
       units = "cm")

load("figures/remake_plots.Rda")
remake_plots[["fig2c"]] =   
  data.frame(
     BAG_corr_mean = df.Test.long$BAG_corr_mean, 
     BAG_corr_change = df.Test.long$BAG_corr_change) 

save(remake_plots, 
     file = "figures/remake_plots.Rda")
  
```


# 4) PCA/feature change
### a) feature change
```{r PCA_ICA - compute feature change}
# open features of interest (if not loaded)
load(file.path(data_folder,"T1w_vars.Rda"))

# open harmonized features
df.Harmonize = read.csv("data/Harmonize.csv", stringsAsFactors = F) %>% 
  mutate(feature = ConensusName,
         modality = as.factor(Stats_file)) %>% 
  dplyr::select(feature, modality)

# compute features change
if (!file.exists(file.path(data_folder,"feature_change.Rda"))) {
  tmp = df.Test %>% 
    dplyr::select(CrossProject_ID,
                  Time,
                   age,
                   BAG_corr_gam,
                   T1w_vars) %>% 
    pivot_longer(T1w_vars, 
                 names_to = "feature", 
                 values_to = "value") %>% 
    group_by(CrossProject_ID, feature) %>% 
    mutate(Time = Time - mean(Time)) %>% 
    do(fit = lm(value ~ Time, data = .))
  
  change = lapply(tmp$fit, function(x) x$coefficients[["Time"]]) %>% simplify2array()
  
  df.slope.feature = data.frame(CrossProject_ID = tmp$CrossProject_ID,
             feature = tmp$feature,
             change = change)
  rm("tmp")
  save(df.slope.feature, file = file.path(data_folder,"feature_change.Rda"))
} else {
  load(file.path(data_folder,"feature_change.Rda"))
}
 

df.feature.change = pivot_wider(df.slope.feature,
                names_from =feature,
                values_from = change)
df.Test.long.f = left_join(df.Test.long, df.feature.change)

tmp = df.Test.long.f %>% dplyr::select(T1w_vars) %>% lapply(., t.test)

df.feature.change = data.frame(feature = names(tmp), 
           p.value = lapply(tmp, function(x) x$p.value) %>% simplify2array(),
           statistic = lapply(tmp, function(x) x$statistic) %>% simplify2array(),
           estimate = lapply(tmp, function(x) x$estimate) %>% simplify2array()) %>% 
  mutate(logp = -log10(p.value),
         logp.sig = if_else(estimate > 0, logp, -logp))

#fdr.thr = -log10(BY(df.feature.change$p.value, alpha = 0.05)$FDR)
bf.thr = -log10(0.05/length(T1w_vars))
sum(df.feature.change$logp > bf.thr)/length(T1w_vars)


df.feature.change = left_join(df.feature.change,df.Harmonize) 

tmp =df.feature.change %>% 
  group_by(modality) %>% 
  summarise(logp =mean(logp)) %>% 
  mutate(order = -rank(logp)) %>% 
  dplyr::select(-logp)
  

df.feature.change = left_join(df.feature.change, tmp) %>% 
arrange(order, -logp) %>% 
 mutate(order2 = 1:length(T1w_vars))

df.feature.change$modality = 
  plyr::mapvalues(df.feature.change$modality,
                from =c("Area_Aparc","GWC_Aparc","Intensity_Aseg","Thickness_Aparc","Volume_Aseg","Volume_Aparc"),
                to = c("area (c)", "gwc (c)","intensity (s)","thickness (c)", "volume (s)", "volume (c)"))


colorscale =c('#1b9e77','#d95f02','#7570b3','#e7298a','#66a61e','#e6ab02')
gs = ggplot(df.feature.change, aes(x = order2, y = logp, group = modality, fill = modality)) +
  geom_point(shape = 21, size = 3, alpha = .1) +
  geom_point(data = df.feature.change %>% filter(logp > bf.thr), shape = 21, size = 3) +
  geom_hline(yintercept = bf.thr, linetype = 3, color ="grey40",size = 1.5) + 
    theme_classic() +
    theme(legend.position = 'none',
        axis.text = element_text(size = 16),
        axis.title = element_text(size = 20),
        axis.title.y = element_text(size = 12),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        legend.title = element_blank(),
        legend.text = element_text(size = 16),
        legend.key.size = unit(5,"point"),
        plot.title = element_text(size = 20, hjust = 0.5, vjust = 0)) +
  xlab("Feature") +
  ylab(expression("Longitudinal Change [-log"[10]*"(p)")) +
  scale_fill_manual(values = colorscale) + 
  scale_y_continuous(trans = "sqrt", breaks = c(1.3,6,20,50,100)) +
  ggtitle("Lifebrain")


ggsave("figures/feature_change_Lifebrain.png", 
       dpi = 500, 
       plot = gs,
       width = 10,
       height = 10, 
       units = "cm")

write.csv(df.feature.change, file = "figures/feature_change_Lifebrain.csv")

df.Test.long %<>% left_join(.,df.Test.long.f %>% dplyr::select(CrossProject_ID, T1w_vars))
```

### b) feature change v. delta_0 and delta_0 change. 

```{r feature and delta}

mod = df.Test.long %>% 
  pivot_longer(T1w_vars, 
               names_to = "feature", 
               values_to = "change") %>% 
  group_by(feature) %>% 
  do(fit = lmer(change ~  BAG_corr_mean+ AgeBsl + sex + scale(eICV) + (1 | cohort), data = ., REML = T))

# UVE - unique variance explained 
mod_uve = df.Test.long %>%
  pivot_longer(T1w_vars,
               names_to = "feature",
               values_to = "change") %>%
  group_by(feature) %>%
  do(fit = lmer(change ~ AgeBsl + sex + scale(eICV) + (1 | cohort), data = ., REML = T))

uve = lapply(mod$fit, function(x) r.squaredGLMM(x)[[1]]) %>% simplify2array() -
      lapply(mod_uve$fit, function(x) r.squaredGLMM(x)[[1]]) %>% simplify2array()
  

df.feature2delta_mean = data.frame(feature = mod$feature, 
           estimate.xgb =lapply(mod$fit, function(x) summary(x)$coefficients["BAG_corr_mean","Estimate"]) %>% simplify2array(),
           p.value.xgb =lapply(mod$fit, function(x) summary(x)$coefficients["BAG_corr_mean","Pr(>|t|)"]) %>% simplify2array(),
            uve.xgb = uve) %>% 
  mutate(logp.xgb = -log10(p.value.xgb),
         logp.sig.xgb = if_else(estimate.xgb > 0, logp.xgb, -logp.xgb))


bf.thr = -log10(0.05/length(T1w_vars))
sum(df.feature2delta_mean$logp.xgb > bf.thr)

df.feature2delta_mean = left_join(df.feature.change %>% dplyr::select(feature,
                                                                      modality,
                                                                      statistic,
                                                                      estimate,
                                                                      logp), 
                                  df.feature2delta_mean, 
                                  by = "feature")

## prepare plot XGB
tmp =df.feature2delta_mean %>% 
  group_by(modality) %>% 
  summarise(logp.xgb =mean(logp.xgb)) %>% 
  mutate(order = -rank(logp.xgb)) %>% 
  dplyr::select(-logp.xgb)
  
df.feature2delta_mean_plot = left_join(df.feature2delta_mean, 
                                  tmp) %>%
  arrange(order, -logp.sig.xgb) %>% 
  mutate(order2 = 1:length(T1w_vars))

colorscale =c('#1b9e77','#d95f02','#7570b3','#e7298a','#66a61e','#e6ab02')
gs = ggplot(df.feature2delta_mean_plot, 
            aes(x = order2, 
                y = logp.sig.xgb, 
                group = modality, 
                fill = modality)) +
  geom_point(shape = 21, 
             size = 3, 
             alpha = .1) +
  geom_point(data = df.feature2delta_mean_plot %>% filter(logp.xgb > bf.thr), 
             shape = 21, 
             size = 3) +
  geom_hline(yintercept = bf.thr, 
             linetype = 3, 
             color ="grey40",
             size = 1.5) + 
  geom_hline(yintercept = -bf.thr, linetype = 3, color ="grey40",size = 1.5) + 
  geom_hline(yintercept = 1.3, linetype = 1, color ="grey40",size = 1.5, alpha =.4) + 
  geom_hline(yintercept = -1.3, linetype = 1, color ="grey40",size = 1.5, alpha =.4) + 
    theme_classic() +
    theme(legend.position = 'none',
        axis.text = element_text(size = 16),
        axis.title = element_text(size = 20),
        axis.title.y = element_text(size = 12),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        legend.title = element_blank(),
        legend.text = element_text(size = 16),
        legend.key.size = unit(5,"point"),
        plot.title = element_text(size = 20, hjust = 0.5, vjust = 0)) + 
  xlab("Feature") +
  ylab(expression("feature change vs. delta"[cross]*" [-log"[10]*"(p)")) +
  scale_fill_manual(values = colorscale) +
  ggtitle("Lifebrain - XGB")

ggsave("figures/feature_change_v_delta_Lifebrain.png", 
       dpi = 500, 
       plot = gs,
       width = 10,
       height = 10, 
       units = "cm")

write.csv(df.feature2delta_mean, file = "figures/feature_change2delta_mean_Lifebrain.csv")
```

### c) feature change v. delta change
```{r feature and delta change}

mod = df.Test.long %>% 
  pivot_longer(T1w_vars, 
               names_to = "feature", 
               values_to = "change") %>% 
  group_by(feature) %>% 
  do(fit = lmer(change ~  BAG_corr_change+ AgeBsl + sex + scale(eICV) + (1 | cohort), data = ., REML = T))

# UVE - unique variance explained 
mod_uve = df.Test.long %>%
  pivot_longer(T1w_vars,
               names_to = "feature",
               values_to = "change") %>%
  group_by(feature) %>%
  do(fit = lmer(change ~ AgeBsl + sex + scale(eICV) + (1 | cohort), data = ., REML = T))

uve = lapply(mod$fit, function(x) r.squaredGLMM(x)[[1]]) %>% simplify2array() -
      lapply(mod_uve$fit, function(x) r.squaredGLMM(x)[[1]]) %>% simplify2array()
  
 
df.feature2delta_change = data.frame(feature = mod$feature, 
           estimate.xgb =lapply(mod$fit, function(x) summary(x)$coefficients["BAG_corr_change","Estimate"]) %>% simplify2array(),
           p.value.xgb =lapply(mod$fit, function(x) summary(x)$coefficients["BAG_corr_change","Pr(>|t|)"]) %>% simplify2array(),
            uve.xgb = uve) %>% 
  mutate(logp.xgb = -log10(p.value.xgb),
         logp.sig.xgb = if_else(estimate.xgb > 0, logp.xgb, -logp.xgb))


bf.thr = -log10(0.05/length(T1w_vars))
sum(df.feature2delta_change$logp.xgb > bf.thr)

df.feature2delta_change = left_join(df.feature.change %>% dplyr::select(feature,
                                                                      modality,
                                                                      statistic,
                                                                      estimate,
                                                                      logp), 
                                  df.feature2delta_change, 
                                  by = "feature")

## prepare plot XGB
tmp =df.feature2delta_change %>% 
  group_by(modality) %>% 
  summarise(logp.xgb =mean(logp.xgb)) %>% 
  mutate(order = -rank(logp.xgb)) %>% 
  dplyr::select(-logp.xgb)
  
df.feature2delta_change_plot = left_join(df.feature2delta_change, 
                                  tmp) %>%
  arrange(order, -logp.sig.xgb) %>% 
  mutate(order2 = 1:length(T1w_vars))

colorscale =c('#1b9e77','#d95f02','#7570b3','#e7298a','#66a61e','#e6ab02')
gs = ggplot(df.feature2delta_change_plot, 
            aes(x = order2, 
                y = logp.sig.xgb, 
                group = modality, 
                fill = modality)) +
  geom_point(shape = 21, 
             size = 3, 
             alpha = .1) +
  geom_point(data = df.feature2delta_change_plot %>% filter(logp.xgb > bf.thr), 
             shape = 21, 
             size = 3) +
  geom_hline(yintercept = bf.thr, 
             linetype = 3, 
             color ="grey40",
             size = 1.5) + 
  geom_hline(yintercept = -bf.thr, linetype = 3, color ="grey40",size = 1.5) + 
  geom_hline(yintercept = 1.3, linetype = 1, color ="grey40",size = 1.5, alpha =.4) + 
  geom_hline(yintercept = -1.3, linetype = 1, color ="grey40",size = 1.5, alpha =.4) + 
    theme_classic() +
    theme(legend.position = 'none',
        axis.text = element_text(size = 16),
        axis.title = element_text(size = 20),
        axis.title.y = element_text(size = 12),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        legend.title = element_blank(),
        legend.text = element_text(size = 16),
        legend.key.size = unit(5,"point"),
        plot.title = element_text(size = 20, hjust = 0.5, vjust = 0)) + 
  xlab("Feature") +
  ylab(expression("feature change vs. delta"[cross]*" [-log"[10]*"(p)")) +
  scale_fill_manual(values = colorscale) +
  ggtitle("Lifebrain - XGB")

ggsave("figures/feature_change_v_delta_change_Lifebrain.png", 
       dpi = 500, 
       plot = gs,
       width = 10,
       height = 10, 
       units = "cm")

write.csv(df.feature2delta_change, file = "figures/feature_change2delta_change_Lifebrain.csv")
```

### d) PCA
```{r PCA on change}
# select only variables with significant change over time
T1w_change = df.feature.change$feature[df.feature.change$logp > bf.thr]

# prcomp
mod.pca = prcomp(df.Test.long[T1w_change], center =F)


df.PCA = data.frame(CrossProject_ID = df.Test.long$CrossProject_ID, 
                    PC1 = mod.pca$x[,1]) %>% 
  mutate(PC1 = if_else(abs(PC1) > 10, NaN, PC1))
summary(mod.pca)$importance[,1:10]

#weights
df.pca.weights = data.frame(
  features = rownames(mod.pca$rotation),
  weight = mod.pca$rotation[,1])
write.csv(df.pca.weights, file = "figures/pca_weights_Lifebrain.csv")

df.Test.long %<>% left_join(., df.PCA)

# basis model
lm.uve = lmer(PC1 ~AgeBsl + sex + scale(eICV)  + (1 | cohort), data = df.Test.long, REML = T)
summary(lm.uve)

# mean xgb
lm.pred = lmer(PC1 ~ BAG_corr_mean + AgeBsl + sex + scale(eICV)  + (1 | cohort), data = df.Test.long, REML = T) 
summary(lm.pred)
r.squaredGLMM(lm.pred) -r.squaredGLMM(lm.uve)

# change xgb
lm.predCh = lmer(PC1 ~ BAG_corr_change + AgeBsl + sex + scale(eICV)  + (1 | cohort), data = df.Test.long, REML = T) 
summary(lm.predCh)
r.squaredGLMM(lm.predCh) -r.squaredGLMM(lm.uve)


```
#### PCA figures
```{r Plot Figure 2a-b. Brain Age to Brain Aging}

#ORANGE SCALE
#yellow scale

#Lifebrain - Extreme Bossting
gs = ggplot(data = df.Test.long,
       mapping = aes(x = BAG_corr_mean,
                     y = PC1)) + 
  geom_point(shape = 21, color = "black", fill = "#EF820D", size = 4, stroke = 1.5) + 
  geom_smooth(method = "lm", 
              color = "#F05E23",
              size = 3) +
  theme_classic() + 
  theme(legend.position = 'none',
        axis.text = element_text(size = 16),
        axis.title = element_text(size = 20),
        #axis.title.y = element_text(size = 12),
        plot.title = element_text(size = 20, hjust = 0.5, vjust = 0)) +
  ylab(expression("Feature change (PC1)")) +
  xlab(expression("Brain Age Delta"[cross])) +
  ggtitle("Lifebrain - XGB")

ggsave("figures/PC1_BrainAgeCross_Lifebrain.png", 
       dpi = 500, 
       plot = gs,
       width = 10,
       height = 10, 
       units = "cm")

# Lifebrain - Extreme Bossting
gs = ggplot(data = df.Test.long,
       mapping = aes(x = BAG_corr_change,
                     y = PC1)) + 
  geom_point(shape = 21, color = "black", fill = "#FCF4A3", size = 4, stroke = 1.5) + 
  geom_smooth(method = "lm", 
              color = "#F8DE7E",
              size = 3) +
  theme_classic() + 
  theme(legend.position = 'none',
        axis.text = element_text(size = 16),
        axis.title.y = element_text(size = 20),
        axis.title = element_text(size = 20),
        plot.title = element_text(size = 20, hjust = 0.5, vjust = 0)) +
  ylab(expression("Feature change (PC1)")) +
  xlab(expression("Brain Age Delta"[long])) +
  ggtitle("Lifebrain - XGB")


ggsave("figures/PC1_BrainAgeLong_Lifebrain.png", 
       dpi = 500, 
       plot = gs,
       width = 10,
       height = 10, 
       units = "cm")

```




### PCA
```{r PCA on change}

#df.Test.long %<>% left_join(.,df.Test.long.f %>% dplyr::select(eid, T1w_vars))
T1w_change = df.feature.change$feature[df.feature.change$logp > bf.thr]

# prcomp
mod.pca = prcomp(df.Test.long[T1w_change], center =F)


data.frame(mod.pca$rotation[,1]) %>% View() # positive values less decline. note different dirsctionality with delta
df.PCA = data.frame(CrossProject_ID = df.Test.long$CrossProject_ID, mod.pca$x[,1:4])
# pc1 = .20, pc2 = .17, pc3 = .04, pc4 = .03

df.Test.long.test = df.Test.long
df.Test.long.test %<>% left_join(., df.PCA)

## uve
lm.uve = lmer(PC1 ~ AgeBsl + sex + scale(eICV)  + (1 | cohort), data = df.Test.long.test, REML = T)

lm.pred = lmer(PC1 ~ BAG_corr_mean + AgeBsl + sex + scale(eICV)  + (1 | cohort), data = df.Test.long.test, REML = T)
summary(lm.pred)
r.squaredGLMM(lm.pred) -r.squaredGLMM(lm.uve)

lm.pred = lmer(PC1 ~ BAG_corr_change + AgeBsl + sex + scale(eICV) + (1 | cohort), data = df.Test.long.test, REML = T)
summary(lm.pred)






```