06_SNR_Results.Rmd

---
title: "SNR Paper Figures"
date: "2024-01-29"
geometry: margin = 1in
output:
  pdf_document: 
    toc: yes
    toc_depth: 4
    fig_caption: yes
  latex_engine: xelatex
---


```{r load-libraries, message=FALSE, echo = FALSE}
knitr::opts_chunk$set(echo = TRUE, warning=FALSE, dev ='pdf', dpi=150, fig.height=4, fig.width = 6, fig.path='/Volumes/Hera/Projects/7TBrainMech/scripts/eeg/Shane/Cortical_SNR_Development/results/rStudioFigsResubmit/')
```

```{r, echo= FALSE, message = FALSE, warning = FALSE}

library(LNCDR)
library(data.table)
library(dplyr)
library(factoextra)
library(ggplot2)
library(e1071)
library(caret)
attach(mtcars)
library(grid)
library(gridExtra)
library(plotrix)
library(mgcv)
library(readxl)
library(lme4)
library(lubridate)
library(checkmate)
library(lmerTest)
library(tidyr)
library(jtools)
library(eegUtils)
library(tvem)
library(interactions)
library(akima)
library(mice)
library(scales)
library(asbio)
library(r2glmm)

tag_facet <- function(p, open = "(", close = ")", tag_pool = letters, x = -Inf, y = Inf, 
                      hjust = -0.5, vjust = 1.5, fontface = 2, ...) {

  gb <- ggplot_build(p)
  lay <- gb$layout$layout
  tags <- cbind(lay, label = paste0(open, tag_pool[lay$PANEL], close), x = x, y = y)
  p + geom_text(data = tags, aes_string(x = "x", y = "y", label = "label"), ..., hjust = hjust, 
                vjust = vjust, fontface = fontface, inherit.aes = FALSE) 
}

```

# SNR across freqs 
```{r SNRacrossFreqs, echo= FALSE, message = FALSE, warning = FALSE, fig.height=10, fig.width=15}

merge7t <- read.csv('/Volumes/Hera/Projects/7TBrainMech/scripts/txt/merged_7t.csv')
sexInfo <- merge7t[c("lunaid","visitno","sex")]
names(sexInfo)[names(sexInfo) == "lunaid"] <- "lunaID"


SNRallChans_outlier <- read.csv('/Volumes/Hera/Projects/7TBrainMech/scripts/eeg/Shane/Cortical_SNR_Development/results/allSubjectsSNR_allChans_allfreqs.csv') %>% merge(., sexInfo, by = c("lunaID", "visitno"))

# pivot to long format
SNRallChans_outlier_long <- pivot_longer(SNRallChans_outlier, cols = c("Total","Evoked", "SNR", "Induced"), names_to = 'measure', values_to = 'value')

# summarize by age group & measurement
SNRallChans_summary_byAge <- SNRallChans_outlier_long %>% 
  mutate(ageGroup = as.factor(cut(age, c(0,17,23,Inf), labels = c('10-16','17-22','23-30')))) %>%
  group_by(freqs, measure, ageGroup, hertz) %>%
  summarize(y = mean(value, na.rm=T),
            sd = sd(value, na.rm=T),
            n = n(),
            se = sd / sqrt(n))

# summarize by age group & measurement
SNRallChans_summary <- SNRallChans_outlier_long %>% 
  mutate(ageGroup = as.factor(cut(age, c(0,17,23,Inf), labels = c('10-16','17-22','23-30')))) %>%
  group_by(freqs, measure, hertz) %>%
  summarize(y = mean(value, na.rm=T),
            sd = sd(value, na.rm=T),
            n = n(),
            se = sd / sqrt(n))


# evoked vs induced
SNRallChans_summary$hertz <- as.factor(SNRallChans_summary$hertz)
levels(SNRallChans_summary$hertz) <- c("20 Hz Clicks", "30 Hz Clicks", "40 Hz Clicks")


lunaize(ggplot(data = SNRallChans_summary %>% filter(freqs <= 50, measure %in% c('Evoked')) %>% filter(hertz == "40 Hz Clicks"), 
       aes(x = freqs, y = y, color = interaction(measure), 
           fill = interaction(measure))) +
  geom_line() +
  geom_ribbon(aes(ymin = y - se, ymax = y + se), alpha = 0.8, linewidth = 1.5) +
  scale_fill_manual(values = c("#E3B23C")) +
  scale_color_manual(values = c("#E3B23C")) +
  scale_y_log10() + 
  theme_bw()) + xlab("Frequences (Hz)") + ylab("Power (\u00b5V\u00b2)") +
   theme(strip.background = element_blank(), strip.placement = "outside", strip.text = element_text(size = 25, vjust = 1.5), text = element_text(size = 40), plot.margin = margin(0.5,0.5,0.5,0.5,"cm"), aspect.ratio = 1, axis.title.y = element_text(vjust = 2), axis.title.x = element_text(vjust = -2), legend.position = "none", axis.text.y = element_blank(), panel.spacing = unit(3, "lines"), axis.ticks.y = element_blank()) 



lunaize(ggplot(data = SNRallChans_summary %>% filter(freqs <= 50, measure %in% c('Induced')) %>% filter(hertz == "40 Hz Clicks"), 
       aes(x = freqs, y = y, color = interaction(measure), 
           fill = interaction(measure))) +
  geom_line() +
  geom_ribbon(aes(ymin = y - se, ymax = y + se), alpha = 0.8, linewidth = 1.5) +
  scale_fill_manual(values = c("plum4")) +
  scale_color_manual(values = c("plum4")) +
  scale_y_log10() + 
  theme_bw()) + xlab("Frequences (Hz)") + ylab("Power (\u00b5V\u00b2)") +
   theme(strip.background = element_blank(), strip.placement = "outside", strip.text = element_text(size = 25, vjust = 1.5), text = element_text(size = 40), plot.margin = margin(0.5,0.5,0.5,0.5,"cm"), aspect.ratio = 1, axis.title.y = element_text(vjust = 2), axis.title.x = element_text(vjust = -2), legend.position = "none", axis.text.y = element_blank(), panel.spacing = unit(3, "lines"), axis.ticks.y = element_blank()) 

lunaize(ggplot(data = SNRallChans_summary %>% filter(freqs <= 50, measure %in% c('SNR')) %>% filter(hertz == "40 Hz Clicks"), 
       aes(x = freqs, y = y, color = interaction(measure), 
           fill = interaction(measure))) +
  geom_line() +
  geom_ribbon(aes(ymin = y - se, ymax = y + se), alpha = 0.8, linewidth = 1.5) +
  scale_fill_manual(values = c("#2E041F")) +
  scale_color_manual(values = c("#2E041F")) +
  scale_y_log10() + 
  theme_bw()) + xlab("Frequencies (Hz)") + ylab("SNR") +
   theme(strip.background = element_blank(), strip.placement = "outside", strip.text.x = element_text(color = "white"), strip.text = element_text(size = 40, vjust = -4), text = element_text(size = 25), plot.margin = margin(0.5,0.5,0.5,0.5,"cm"), aspect.ratio = 1, axis.title.y = element_text(vjust = 2), axis.title.x = element_text(vjust = -2), legend.position = "none", axis.text.y = element_blank(), panel.spacing = unit(3, "lines"), axis.ticks.y = element_blank()) 

```

# Topo plots
```{r topo_plots, echo= FALSE, message = FALSE, warning = FALSE, fig.height=10, fig.width=25}

SNRchanLocs <- read.csv('/Volumes/Hera/Projects/7TBrainMech/scripts/eeg/Shane/Cortical_SNR_Development/results/SNRchanlocs.csv')

evokedchanLocs <- read.csv('/Volumes/Hera/Projects/7TBrainMech/scripts/eeg/Shane/Cortical_SNR_Development/results/evokedchanLocs.csv')
evokedchanLocs$measure <- "Evoked"

inducedchanLocs <- read.csv('/Volumes/Hera/Projects/7TBrainMech/scripts/eeg/Shane/Cortical_SNR_Development/results/inducedchanLocs.csv') %>% select(-measure)
inducedchanLocs$measure <- "Spontaneous"

allmeasuresChanLocs <- rbind(SNRchanLocs, evokedchanLocs) %>% rbind(., inducedchanLocs)

lunaize(ggplot(allmeasuresChanLocs %>% filter(measure == "Spontaneous"), aes(x = -Y, y = X, fill = PC1, z = PC1, label = labels)) + 
          geom_topo(chan_markers = "text", interpolate = F, interp_limit = "head") +
          scale_fill_gradient2(low="white", mid="white", high="plum4")) +
   theme(strip.background = element_blank(), strip.placement = "outside", strip.text = element_text(size = 40), text = element_text(size = 25), plot.margin = margin(0.5,0.5,0.5,0.5,"cm"), aspect.ratio = 1, axis.title.y = element_text(vjust = 2), axis.title.x = element_text(vjust = -2), panel.spacing = unit(3, "lines"))

lunaize(ggplot(allmeasuresChanLocs %>% filter(measure == "Evoked"), aes(x = -Y, y = X, fill = PC1, z = PC1, label = labels)) + 
          geom_topo(chan_markers = "text", interpolate = F, interp_limit = "head") +
          scale_fill_gradient2(low="white", mid="white", high="#E3B23C")) +
   theme(strip.background = element_blank(), strip.placement = "outside", strip.text = element_text(size = 40), text = element_text(size = 25), plot.margin = margin(0.5,0.5,0.5,0.5,"cm"), aspect.ratio = 1, axis.title.y = element_text(vjust = 2), axis.title.x = element_text(vjust = -2), panel.spacing = unit(3, "lines"))



lunaize(ggplot(allmeasuresChanLocs %>% filter(measure == "SNR"), aes(x = -Y, y = X, fill = PC1, z = PC1, label = labels)) + 
          geom_topo(chan_markers = "text", interpolate = F, interp_limit = "head") +
          scale_fill_gradient2(low="white", mid="white", high="#2E041F")) +
   theme(strip.background = element_blank(), strip.placement = "outside", strip.text = element_text(size = 40), text = element_text(size = 25), plot.margin = margin(0.5,0.5,0.5,0.5,"cm"), aspect.ratio = 1, axis.title.y = element_text(vjust = 2), axis.title.x = element_text(vjust = -2), panel.spacing = unit(3, "lines"))

```

# Age plots
```{r age_plots, echo= FALSE, message = FALSE, warning = FALSE, fig.height=18, fig.width=15}

allSNRmeasures <- read.csv('/Volumes/Hera/Projects/7TBrainMech/scripts/eeg/Shane/Cortical_SNR_Development/results/SNRmeasures_PCAvalues.csv')
allSNRmeasures$Stim <- "40 Hz"
allSNRmeasures$lunaID[allSNRmeasures$lunaID == "11748"] <- "11515"
allSNRmeasures$lunaID[allSNRmeasures$lunaID == "11665"] <- "11390"
allSNRmeasures <- merge(allSNRmeasures, sexInfo, by = c("lunaID", "visitno"), all = T)


allSNRmeasures$sex[allSNRmeasures$lunaID == "11653" & is.na(allSNRmeasures$sex)] <- "F"
allSNRmeasures$sex[allSNRmeasures$lunaID == "11690" & is.na(allSNRmeasures$sex)] <- "F"
allSNRmeasures$sex[allSNRmeasures$lunaID == "11812" & is.na(allSNRmeasures$sex)] <- "F"



allSNRmeasures$sex[allSNRmeasures$lunaID == "11515" & is.na(allSNRmeasures$sex)] <- "F"
allSNRmeasures$sexFac <- factor(allSNRmeasures$sex, levels = c("M", "F"), ordered = T)


#allSNRmeasures$inverseAge <- 1/(allSNRmeasures$age)

#SNR 
lunaize(ggplot(data = allSNRmeasures %>% filter(measure == "SNR"), aes(x = age, y = pc1)) +
          geom_point(size = 5,alpha=.8, color = "#2E041F") +
            geom_line(aes(group= interaction(lunaID)), alpha = 0.1) + 
            geom_smooth(method=mgcv::"gam", formula = y ~ s(x, k = 3, fx = T), alpha=0.4, size=2, color = "#2E041F")) + 
  xlab("Age") + ylab("SNR PC1") +
  geom_text(aes(y = -0.9, x = 34.5, label = "*"), size = 15, color="#2E041F") + 
   theme(text = element_text(size = 40), plot.margin = margin(0.5,0.5,0.5,0.5,"cm"), aspect.ratio = 1,
        axis.title.y = element_text(vjust = 2), 
        axis.title.x = element_text(vjust = -2)) 
 
gam.model <-  mgcv::gamm(pc1 ~ s(age, k = 3) + sexFac, data = allSNRmeasures %>% filter(measure == "SNR"), random=list(lunaID=~1))
summary(gam.model$gam)

# sex interactions
gam.model <-  mgcv::gamm(pc1 ~ s(age, by = sexFac) + sexFac, data = allSNRmeasures %>% filter(measure == "SNR"), random=list(lunaID=~1))
summary(gam.model$gam)
 
#evoked
lunaize(ggplot(data = allSNRmeasures %>% filter(measure == "Evoked"), aes(x = age, y = pc1)) +
          geom_point(size =5,alpha=.8, color = "#E3B23C") +
            geom_line(aes(group= interaction(lunaID)), alpha = 0.1) + 
            geom_smooth(method=mgcv::"gam", formula = y ~ s(x, k = 3, fx = T), alpha=0.4, size=2, color = "#E3B23C")) + 
  xlab("Age") + ylab("Evoked PC1") + 
   theme(text = element_text(size = 40), plot.margin = margin(0.5,0.5,0.5,0.5,"cm"), aspect.ratio = 1,
        axis.title.y = element_text(vjust = 2), 
        axis.title.x = element_text(vjust = -2)) 

gam.model <-  mgcv::gamm(pc1 ~ s(age, k = 3) + sexFac, data = allSNRmeasures %>% filter(measure == "Evoked"), random=list(lunaID=~1))
summary(gam.model$gam)

# sex interactions
gam.model <-  mgcv::gamm(pc1 ~ s(age, by = sexFac) + sexFac, data = allSNRmeasures %>% filter(measure == "Evoked"), random=list(lunaID=~1))
summary(gam.model$gam)


#induced
lunaize(ggplot(data = allSNRmeasures %>% filter(measure == "Induced"), aes(x = age, y = pc1)) +
          geom_point(size =5,alpha=.8, color = "plum4") +
            geom_line(aes(group= interaction(lunaID)), alpha = 0.1) + 
            geom_smooth(method=mgcv::"gam", formula = y ~ s(x, k = 3, fx = T), alpha=0.4, size=2, color="plum4")) + 
  xlab("Age") + ylab("Spontaneous PC1") +
  geom_text(aes(y = -3, x = 35.5, label = "***"), size = 15, color="plum4") + 
   theme(text = element_text(size = 40), plot.margin = margin(0.5,0.5,0.5,0.5,"cm"), aspect.ratio = 1,
        axis.title.y = element_text(vjust = 2), 
        axis.title.x = element_text(vjust = -2)) 
  
gam.model <-  mgcv::gamm(pc1 ~ s(age, k = 3) + sexFac, data = allSNRmeasures %>% filter(measure == "Induced"), random=list(lunaID=~1))
summary(gam.model$gam)

# sex interactions
gam.model <-  mgcv::gamm(pc1 ~ s(age, by = sexFac) + sexFac, data = allSNRmeasures %>% filter(measure == "Induced"), random=list(lunaID=~1))
summary(gam.model$gam)

lunaize(ggplot(data = allSNRmeasures %>% filter(measure == "Induced"), aes(x = age, y = pc1)) +
          geom_point(size =5,alpha=.8, color = "plum4") +
            geom_line(aes(group= interaction(lunaID)), alpha = 0.1) + 
            geom_smooth(aes(group=sex, color = sex),method=mgcv::"gam", formula = y ~ s(x, k = 3, fx = T), alpha=0.4, size=2)) + 
  xlab("Age") + ylab("Spontaneous PC1") +
   theme(text = element_text(size = 40), plot.margin = margin(0.5,0.5,0.5,0.5,"cm"), aspect.ratio = 1,
        axis.title.y = element_text(vjust = 2), 
        axis.title.x = element_text(vjust = -2)) 

```


# Load MRSI 
```{r, echo= FALSE, message = FALSE, warning = FALSE, fig.height=10, fig.width=10}

MRS <- read.csv('/Volumes/Hera/Projects/7TBrainMech/scripts/eeg/Shane/Cortical_SNR_Development/results/allSubjectsDLPFCMRSMeasures_20240223.csv') %>% select(-X)
colnames(MRS) <- c("lunaID","visitno","age","region","GABA","Glu", "GluGABAimbalance", "GluGABAimbalanceABS", "Ratio")

# merge PCA components and MRS
allmeasures_mrs <- merge(allSNRmeasures, MRS, by = c("lunaID", "visitno", "age"))
allmeasures_mrs$inverseAge <- 1/allmeasures_mrs$age
```

### Glu
```{r Glu , echo= FALSE, message = FALSE, warning = FALSE, fig.height=10, fig.width=10}

#SNR 

lunaize(ggplot(data = allmeasures_mrs %>% filter(measure == 'SNR'), aes(x = Glu, y = pc1)) +
            geom_point(size = 3,alpha=1, color = "#2E041F") +
          geom_line(aes(group= interaction(lunaID, region), color = "#2E041F"), alpha = 0.3, size = 2) +
          geom_smooth(aes(group = 1), method="lm", alpha = 0.4, color="black",size=3)) + 
  xlab("Glutamate") + ylab("SNR PC1") + 
         scale_color_manual(values=c("#2E041F")) +
   theme(text = element_text(size = 40), plot.margin = margin(0.5,0.5,0.5,0.5,"cm"), aspect.ratio = 1,
        axis.title.y = element_text(vjust = 2), 
        axis.title.x = element_text(vjust = -2), legend.position = "none") 

lm.model <-  lmer(Glu ~ age + pc1 + sexFac + region + (1|lunaID), data = allmeasures_mrs %>% filter(measure == "SNR"))
summary(lm.model)

#age interaction
lm.model <-  lmer(Glu ~ age *pc1 + region + sexFac + (1|lunaID), data = allmeasures_mrs %>% filter(measure == "SNR"))
summary(lm.model)


#spontaneous
 lunaize(ggplot(data = allmeasures_mrs %>% filter(measure == 'Induced'), aes(x = Glu, y = pc1)) +
           geom_point(size = 3,alpha=1, color = "plum4") +
          geom_line(aes(group= interaction(lunaID, region), color = "plum4"), alpha = 0.3, size = 2) +
          geom_smooth(aes(group = 1), method="lm", alpha = 0.4, color="black",size=3)) + 
  xlab("Glutamate") + ylab("Spontaneous PC1") + 
         scale_color_manual(values=c("plum4")) +
   theme(text = element_text(size = 40), plot.margin = margin(0.5,0.5,0.5,0.5,"cm"), aspect.ratio = 1,
        axis.title.y = element_text(vjust = 2), 
        axis.title.x = element_text(vjust = -2), legend.position = "none") +
    geom_text(aes(y = -1, x = 2.4, label = "*"), size = 15, color="plum4") 


lm.model <-  lmer(Glu ~ age + pc1 + region + sexFac + (1|lunaID), data = allmeasures_mrs %>% filter(measure == "Induced"))
summary(lm.model)

#age interaction
lm.model <-  lmer(Glu ~ age *pc1 + region + sexFac + (1|lunaID), data = allmeasures_mrs %>% filter(measure == "Induced"))
summary(lm.model)


```

### GABA
```{r GABA , echo= FALSE, message = FALSE, warning = FALSE, fig.height=10, fig.width=10}

#SNR 

lunaize(ggplot(data = allmeasures_mrs %>% filter(measure == 'SNR'), aes(x = GABA, y = pc1)) +
            geom_point(size = 2,alpha=0.8, color = "#2E041F") +
          geom_line(aes(group= interaction(lunaID, region), color = "#2E041F"), alpha = 0.2) +
          geom_smooth(aes(group = 1), method="lm", alpha = 0.4, color="black",size=2)) + 
  xlab("GABA") + ylab("SNR PC1") + 
         scale_color_manual(values=c("#2E041F")) +
   theme(text = element_text(size = 40), plot.margin = margin(0.5,0.5,0.5,0.5,"cm"), aspect.ratio = 1,
        axis.title.y = element_text(vjust = 2), 
        axis.title.x = element_text(vjust = -2), legend.position = "none") 

lm.model <-  lmer(GABA ~ age + pc1 + sexFac + region + (1|lunaID), data = allmeasures_mrs %>% filter(measure == "SNR"))
summary(lm.model)

#age interaction
lm.model <-  lmer(GABA ~ age *pc1 +region + sexFac + (1|lunaID), data = allmeasures_mrs %>% filter(measure == "SNR"))
summary(lm.model)



#spontaneous
 lunaize(ggplot(data = allmeasures_mrs %>% filter(measure == 'Induced'), aes(x = GABA, y = pc1)) +
            geom_point(size = 2,alpha=0.8, color = "plum4") +
          geom_line(aes(group= interaction(lunaID, region), color = "plum4"), alpha = 0.2) +
          geom_smooth(aes(group = 1), method="lm", alpha = 0.4, color="black",size=2)) + 
  xlab("GABA") + ylab("Spontaneous PC1") + 
         scale_color_manual(values=c("plum4")) +
   theme(text = element_text(size = 40), plot.margin = margin(0.5,0.5,0.5,0.5,"cm"), aspect.ratio = 1,
        axis.title.y = element_text(vjust = 2), 
        axis.title.x = element_text(vjust = -2), legend.position = "none") 

 
lm.model <-  lmer(GABA ~ age + pc1 + sexFac + region + (1|lunaID), data = allmeasures_mrs %>% filter(measure == "Induced"))
summary(lm.model)

#age interaction
lm.model <-  lmer(GABA ~ age *pc1 +region + (1|lunaID), data = allmeasures_mrs %>% filter(measure == "Induced"))
summary(lm.model)


```



### Imbalance
```{r Imbalance, echo= FALSE, message = FALSE, warning = FALSE, fig.height=10, fig.width=18}

# SNR pc1

 lunaize(ggplot(data = allmeasures_mrs %>% filter(measure == 'SNR'), aes(x = sqrt(GluGABAimbalanceABS), y = pc1)) +
            geom_point(size = 3,alpha=1, color = "#2E041F") +
          geom_line(aes(group= interaction(lunaID, region), color = "#2E041F"), alpha = 0.3, size = 2) +
          geom_smooth(aes(group = 1), method="lm", alpha = 0.4, color="black",size=3)) + 
  xlab("Glu GABA Asymmetry") + ylab("SNR PC1") + 
         scale_color_manual(values=c("#2E041F")) +
   theme(text = element_text(size = 40), plot.margin = margin(0.5,0.5,0.5,0.5,"cm"), aspect.ratio = 1,
        axis.title.y = element_text(vjust = 2), 
        axis.title.x = element_text(vjust = -2), legend.position = "none")+ 
  geom_text(aes(x = 0.98, y = -1.3, label = "*"), size = 15, color="black") 


lm.model <-  lmer(sqrt(GluGABAimbalanceABS) ~ age + pc1 + sexFac + region + (1|lunaID), data = allmeasures_mrs %>% filter(measure == "SNR"))
summary(lm.model)
  
#age interaction
lm.model <-  lmer(pc1 ~ age*sqrt(GluGABAimbalanceABS) + sexFac + region+ (1|lunaID), data = allmeasures_mrs %>% filter(measure == "SNR"))
summary(lm.model)

# spontaneous pc1

 lunaize(ggplot(data = allmeasures_mrs %>% filter(measure == 'Induced'), aes(x = sqrt(GluGABAimbalanceABS), y = pc1)) +
               geom_point(size = 3,alpha=1, color = "plum4") +
          geom_line(aes(group= interaction(lunaID, region), color = "plum4"), alpha = 0.3, size = 2) +
          geom_smooth(aes(group = 1), method="lm", alpha = 0.4, color="black",size=3)) +  
  xlab("Glu GABA Asymmetry") + ylab("Spontaneous PC1") + 
         scale_color_manual(values=c("plum4")) +
   theme(text = element_text(size = 40), plot.margin = margin(0.5,0.5,0.5,0.5,"cm"), aspect.ratio = 1,
        axis.title.y = element_text(vjust = 2), 
        axis.title.x = element_text(vjust = -2), legend.position = "none")


lm.model <-  lmer(sqrt(GluGABAimbalanceABS) ~ age + pc1 + sexFac + region + (1|lunaID), data = allmeasures_mrs %>% filter(measure == "Induced"))
summary(lm.model)
  
#age interaction
lm.model <-  lmer(pc1 ~ age*sqrt(GluGABAimbalanceABS) + sexFac + region+ (1|lunaID), data = allmeasures_mrs %>% filter(measure == "Induced"))
summary(lm.model)




```

### Ratio
```{r ratio, echo= FALSE, message = FALSE, warning = FALSE, fig.height=10, fig.width=18}

# SNR pc1

 lunaize(ggplot(data = allmeasures_mrs %>% filter(measure == 'SNR'), aes(x = sqrt(Ratio), y = pc1)) +
            geom_point(size = 2,alpha=0.8, color = "#2E041F") +
          geom_line(aes(group= interaction(lunaID, region), color = "#2E041F"), alpha = 0.2) +
          geom_smooth(aes(group = 1), method="lm", alpha = 0.4, color="black",size=2)) + 
  xlab("Glu GABA Ratio") + ylab("SNR PC1") + 
         scale_color_manual(values=c("#2E041F")) +
   theme(text = element_text(size = 40), plot.margin = margin(0.5,0.5,0.5,0.5,"cm"), aspect.ratio = 1,
        axis.title.y = element_text(vjust = 2), 
        axis.title.x = element_text(vjust = -2), legend.position = "none")


lm.model <-  lmer(sqrt(Ratio) ~ age + pc1 + region + (1|lunaID), data = allmeasures_mrs %>% filter(measure == "SNR"))
summary(lm.model)
  
#age interaction
lm.model <-  lmer(pc1 ~ age*sqrt(Ratio) + region+ (1|lunaID), data = allmeasures_mrs %>% filter(measure == "SNR"))
summary(lm.model)


# spontaneous pc1

 lunaize(ggplot(data = allmeasures_mrs %>% filter(measure == 'Induced'), aes(x = sqrt(Ratio), y = pc1)) +
            geom_point(size = 2,alpha=0.8, color = "plum4") +
          geom_line(aes(group= interaction(lunaID, region), color = "plum4"), alpha = 0.2) +
          geom_smooth(aes(group = 1), method="lm", alpha = 0.4, color="black",size=2)) + 
  xlab("Glu GABA Ratio") + ylab("Spontaneous PC1") + 
         scale_color_manual(values=c("plum4")) +
   theme(text = element_text(size = 40), plot.margin = margin(0.5,0.5,0.5,0.5,"cm"), aspect.ratio = 1,
        axis.title.y = element_text(vjust = 2), 
        axis.title.x = element_text(vjust = -2), legend.position = "none")

lm.model <-  lmer(sqrt(Ratio) ~ age + pc1 + region + (1|lunaID), data = allmeasures_mrs %>% filter(measure == "Induced"))
summary(lm.model)
  
#age interaction
lm.model <-  lmer(pc1 ~ age*sqrt(Ratio) + region+ (1|lunaID), data = allmeasures_mrs %>% filter(measure == "Induced"))
summary(lm.model)


```


# Behavior
```{r bestError, echo= FALSE, message = FALSE, warning = FALSE, fig.height=15, fig.width=30}
behav <- read.csv("/Users/shanemckeon/Library/CloudStorage/OneDrive-UniversityofPittsburgh/Lab/Projects/FOOOF/R_Markdown/allSubjectsBehavior.csv")
names(behav)[names(behav) == "luna"] <- "lunaID"

SNRbehavior <-  merge(allSNRmeasures, behav, by= c("lunaID", "visitno", "age")) %>% select(-X.x, -X.y)

SNRbehavior$measure <- as.factor(SNRbehavior$measure)
levels(SNRbehavior$measure) <- c("Evoked", "Spontaneous", "SNR")


SNRbehavior$zscore_Acc <- zscore(SNRbehavior$absBestError)
SNRbehavior$zscore_Accsd <- zscore(SNRbehavior$absBestError_sd)
SNRbehavior$zscore_Lat <- zscore(SNRbehavior$mgsLatency)
SNRbehavior$zscore_Latsd <- zscore(SNRbehavior$mgsLatency_sd)


SNRbehavior_long <- SNRbehavior %>%
  pivot_longer(
    cols = starts_with("zscore_"),
    names_to = "behavioral",
    names_prefix = "zscore_",
    values_to = "value"
  )

lunaize(ggplot(data = SNRbehavior_long %>% 
  mutate(ageGroup = as.factor(cut(age, c(0,15,21,Inf), labels = c('10-14','15-20','21-30')))), aes(y = value, x = age, color = behavioral)) +
          geom_smooth(aes(group = behavioral), method="lm",formula = 'y~I(1/x)', alpha = 0.2,size=1)) + 
  ylab("Behavioral Measure Zscore") + xlab("Age") +      
  scale_color_manual(values=c("#E3B23C", "plum4","#2E041F","black"))+
   theme(strip.background = element_blank(), strip.placement = "outside", strip.text = element_text(size = 30), text = element_text(size = 40), plot.margin = margin(0.5,0.5,0.5,0.5,"cm"),  legend.position = "right", aspect.ratio = 1, axis.title.y = element_text(vjust = 2),
        axis.title.x = element_text(vjust = -2), panel.spacing = unit(3, "lines"))


# accuracy 
lunaize(ggplot(data = SNRbehavior %>% 
  mutate(ageGroup = as.factor(cut(age, c(0,15,21,Inf), labels = c('10-14','15-20','21-30')))), aes(y = pc1, x = absBestError, color = measure)) +
          geom_point(size = 2,alpha=0.8) +
          geom_line(aes(group= interaction(lunaID)), alpha = 0.2) +
          geom_smooth(aes(group = 1, color = measure), method="lm", alpha = 0.4,size=2)) + 
  ylab("PC1") + xlab("Accuracy (degs)") + facet_wrap(~measure)+      
  scale_color_manual(values=c("#E3B23C", "plum4","#2E041F"))+
   theme(strip.background = element_blank(), strip.placement = "outside", strip.text = element_text(size = 30), text = element_text(size = 40), plot.margin = margin(0.5,0.5,0.5,0.5,"cm"),  legend.position = "none", aspect.ratio = 1, axis.title.y = element_text(vjust = 2),
        axis.title.x = element_text(vjust = -2), panel.spacing = unit(3, "lines"))

## vs SNR
lm.model <-  lmer(pc1 ~ age+ absBestError + sexFac + (1|lunaID), data = SNRbehavior %>% filter(measure == 'SNR'))
summary(lm.model)

lm.model <-  lmer(pc1 ~ age*absBestError + sexFac + (1|lunaID), data = SNRbehavior %>% filter(measure == 'SNR'))
summary(lm.model)

## vs evoked
lm.model <-  lmer(pc1 ~ age+ absBestError + sexFac + (1|lunaID), data = SNRbehavior %>% filter(measure == 'Evoked'))
summary(lm.model)

lm.model <-  lmer(pc1 ~ age*absBestError + sexFac + (1|lunaID), data = SNRbehavior %>% filter(measure == 'Evoked'))
summary(lm.model)

## vs spontaneous
lm.model <-  lmer(pc1 ~ age+ absBestError + sexFac + (1|lunaID), data = SNRbehavior %>% filter(measure == 'Spontaneous'))
summary(lm.model)

lm.model <-  lmer(pc1 ~ age*absBestError + sexFac + (1|lunaID), data = SNRbehavior %>% filter(measure == 'Spontaneous'))
summary(lm.model)

# accuracy var

lunaize(ggplot(data = SNRbehavior %>% 
  mutate(ageGroup = as.factor(cut(age, c(0,15,21,Inf), labels = c('10-14','15-20','21-30')))), aes(y = pc1, x = absBestError_sd, color = measure)) + geom_point(size = 2,alpha=0.8) +
          geom_line(aes(group= interaction(lunaID)), alpha = 0.2) +
          geom_smooth(aes(group = 1, color = measure), method="lm", alpha = 0.4,size=2)) + 
  ylab("PC1") + xlab("Accuracy Var (degs)") + facet_wrap(~measure)+      
  scale_color_manual(values=c("#E3B23C", "plum4","#2E041F"))+
   theme(strip.background = element_blank(), strip.placement = "outside", strip.text = element_text(size = 30), text = element_text(size = 40), plot.margin = margin(0.5,0.5,0.5,0.5,"cm"),  legend.position = "none", aspect.ratio = 1, axis.title.y = element_text(vjust = 2),
        axis.title.x = element_text(vjust = -2), panel.spacing = unit(3, "lines"))

## vs SNR
lm.model <-  lmer(pc1 ~ age+ absBestError_sd + sexFac + (1|lunaID), data = SNRbehavior %>% filter(measure == 'SNR'))
summary(lm.model)

lm.model <-  lmer(pc1 ~ age*absBestError_sd + sexFac + (1|lunaID), data = SNRbehavior %>% filter(measure == 'SNR'))
summary(lm.model)

## vs evoked
lm.model <-  lmer(pc1 ~ age+ absBestError_sd + sexFac + (1|lunaID), data = SNRbehavior %>% filter(measure == 'Evoked'))
summary(lm.model)

lm.model <-  lmer(pc1 ~ age*absBestError_sd + sexFac + (1|lunaID), data = SNRbehavior %>% filter(measure == 'Evoked'))
summary(lm.model)

## vs spontaneous
lm.model <-  lmer(pc1 ~ age+ absBestError_sd + sexFac + (1|lunaID), data = SNRbehavior %>% filter(measure == 'Spontaneous'))
summary(lm.model)

lm.model <-  lmer(pc1 ~ age*absBestError_sd + sexFac + (1|lunaID), data = SNRbehavior %>% filter(measure == 'Spontaneous'))
summary(lm.model)


lunaize(ggplot(data = SNRbehavior %>% 
  mutate(ageGroup = as.factor(cut(age, c(0,15,21,Inf), labels = c('10-14','15-20','21-30')))), aes(y = pc1, x = mgsLatency, color = measure)) + geom_point(size = 2,alpha=0.8) +
          geom_line(aes(group= interaction(lunaID)), alpha = 0.2) +
          geom_smooth(aes(group = 1, color = measure), method="lm", alpha = 0.4,size=2)) + 
  ylab("PC1") + xlab("MGS Latency (s)") + facet_wrap(~measure)+      
  scale_color_manual(values=c("#E3B23C", "plum4","#2E041F"))+
   theme(strip.background = element_blank(), strip.placement = "outside", strip.text = element_text(size = 30), text = element_text(size = 40), plot.margin = margin(0.5,0.5,0.5,0.5,"cm"),  legend.position = "none", aspect.ratio = 1, axis.title.y = element_text(vjust = 2),
        axis.title.x = element_text(vjust = -2), panel.spacing = unit(3, "lines"))+ scale_x_continuous(labels = label_number(accuracy = 0.01))



lm.model <-  lmer(pc1 ~ age+ mgsLatency_sd + sexFac + (1|lunaID), data = SNRbehavior %>% filter(measure == 'Evoked'))
summary(lm.model)

lm.model <-  lmer(mgsLatency_sd ~ age+ pc1 + sexFac + (1|lunaID), data = SNRbehavior %>% filter(measure == 'Evoked'))
summary(lm.model)

p.adjust(0.023, method = 'bonferroni', n=3)


lm.model <-  lmer(pc1 ~ age*mgsLatency_sd + sexFac + (1|lunaID), data = SNRbehavior %>% filter(measure == 'Evoked'))
summary(lm.model)



lunaize(ggplot(data = SNRbehavior %>% 
  mutate(ageGroup = as.factor(cut(age, c(0,15,21,Inf), labels = c('10-14','15-20','21-30')))), aes(y = pc1, x = mgsLatency_sd, color = measure)) + geom_point(size = 2,alpha=0.8) +
          geom_line(aes(group= interaction(lunaID)), alpha = 0.2) +
          geom_smooth(aes(group = 1, color = measure), method="lm", alpha = 0.4,size=2)) + 
  ylab("PC1") + xlab("MGS Latency Var (s)") + facet_wrap(~measure)+      
  scale_color_manual(values=c("#E3B23C", "plum4","#2E041F"))+
   theme(strip.background = element_blank(), strip.placement = "outside", strip.text = element_text(size = 30), text = element_text(size = 40), plot.margin = margin(0.5,0.5,0.5,0.5,"cm"),  legend.position = "none", aspect.ratio = 1, axis.title.y = element_text(vjust = 2),
        axis.title.x = element_text(vjust = -2), panel.spacing = unit(3, "lines")) + scale_x_continuous(labels = label_number(accuracy = 0.01))



lm.model <-  lmer(pc1 ~ age+ mgsLatency_sd + sexFac + (1|lunaID), data = SNRbehavior %>% filter(measure == 'SNR'))
summary(lm.model)

lm.model <-  lmer(pc1 ~ age*mgsLatency_sd + sexFac + (1|lunaID), data = SNRbehavior %>% filter(measure == 'SNR'))
summary(lm.model)

#spontaneous
lm.model <-  lmer(pc1 ~ age+ mgsLatency_sd + sexFac + (1|lunaID), data = SNRbehavior %>% filter(measure == 'Spontaneous'))
summary(lm.model)

lm.model <-  lmer(pc1 ~ age*mgsLatency_sd + sexFac + (1|lunaID), data = SNRbehavior %>% filter(measure == 'Spontaneous'))
summary(lm.model)



# mgs latency

#evoked
lm.model <-  lmer(pc1 ~ age+ mgsLatency + sexFac + (1|lunaID), data = SNRbehavior %>% filter(measure == 'Evoked'))
summary(lm.model)

lm.model <-  lmer(pc1 ~ age*mgsLatency + sexFac + (1|lunaID), data = SNRbehavior %>% filter(measure == 'Evoked'))
summary(lm.model)


#spontaneous
lm.model <-  lmer(pc1 ~ age+ mgsLatency + sexFac + (1|lunaID), data = SNRbehavior %>% filter(measure == 'Spontaneous'))
summary(lm.model)

lm.model <-  lmer(pc1 ~ age*mgsLatency + sexFac + (1|lunaID), data = SNRbehavior %>% filter(measure == 'Spontaneous'))
summary(lm.model)


#SNR
lm.model <-  lmer(pc1 ~ age+ mgsLatency + sexFac + (1|lunaID), data = SNRbehavior %>% filter(measure == 'SNR'))
summary(lm.model)

lm.model <-  lmer(pc1 ~ age*mgsLatency + sexFac + (1|lunaID), data = SNRbehavior %>% filter(measure == 'SNR'))
summary(lm.model)


#cowplot::plot_grid(bestErrorSD, mgsLat, ncol = 1)

```
# Behavior Loops 
```{r behavLoop, echo= FALSE, message = FALSE, warning = FALSE, fig.height=20, fig.width=22}

df_long <- SNRbehavior %>%
  pivot_longer(
    cols = starts_with("pc"),
    names_to = "pc",
    values_to = "pc_value"
  )

# Pivot wider to spread the measure values
df_wide <- df_long %>%
  unite("pc_measure", pc, measure, sep = "_") %>%
  pivot_wider(
    names_from = pc_measure,
    values_from = pc_value
  )


behavVars <- c('absBestError','mgsLatency_sd')
PCAVars <- c('pc1_Spontaneous','pc1_Evoked', 'pc1_SNR')

# main effect
output <- c()  
for (behavVar in behavVars) {
  for (PCAVar in PCAVars) {
   
    model <- paste0(behavVar, ' ~ ', PCAVar, '+ age + (1|lunaID)')
    
    model.out <- summary(lmerTest::lmer(model, data = df_wide))  
    b <- model.out$coefficients[2,1]
    t <- model.out$coefficients[2,4]
    p <- model.out$coefficients[2,5]
    pcor <- p.adjust((p), method = "bonferroni", n = 3)

    lm.with <- lmerTest::lmer(model, data = df_wide)

    out <- r2beta(lm.with, partial = TRUE, method = "sgv", data = NULL)

    R2 <- out$Rsq[3]
    
    output <- rbind(output, data.frame(PCAVar, behavVar, R2, b, t, pcor))   
  }
}  

```