EEG_bisection/derivates/EEG_analysis.Rmd

---
title: "Electrophysiological signatures of temporal context in the bisection task"
author: "Cemre Baykan, Xiuna Zhu, Artyom Zinchenko, Hermann Muller, Zhuanghua Shi"
date: "February 2023"
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
source('dataana.R')
```

#EXPERIMENT 1

```{r}
# to see GM and AM of the sets
mDist_exp= eeg_dat_all %>% dplyr::group_by(Exp, cond) %>%
  dplyr::summarise(n= n(),
                   am_Dur = mean(targetDur),
                   sd_Dur = sd(targetDur),
                   gm_Dur = gm_mean(targetDur),
                   se_Dur = sd(targetDur)/sqrt(n-1))
mDist_exp
```

```{r}
cnv_fit <- function(df){
  fit = lm(CNV ~ t, data = df)
  return(tibble(B = fit$coefficients[1], slope= fit$coefficients[2]))
}
```

# 1.CNV
## Slopes

```{r}
allAverageDat_ica_exp1_ns_n <- allAverageDat_exp1 %>% filter(cond == 'NS')
allAverageDat_ica_exp1_ns_n$targetDur = allAverageDat_ica_exp1_ns_n$targetDur /1000
allAverageDat_ica_exp1_ns_n$t = allAverageDat_ica_exp1_ns_n$t /1000
allAverageDat_ica_exp1_ns_n %>% filter(t <= 0.65& t >= 0.25) %>% 
  group_by(SubName) %>% nest() %>%
  mutate(model = map(data, cnv_fit)) %>% 
  unnest(model, .drop = TRUE) -> lm_exp1_ns_slope1
lm_exp1_ns_slope1$cond <- "NS"

t.test(lm_exp1_ns_slope1$slope, mu = 0, alternative = "two.sided")
```


```{r}
#in ps
allAverageDat_ica_exp1_ps_n <- allAverageDat_exp1 %>% filter(cond == 'PS')
allAverageDat_ica_exp1_ps_n$targetDur = allAverageDat_ica_exp1_ps_n$targetDur /1000
allAverageDat_ica_exp1_ps_n$t = allAverageDat_ica_exp1_ps_n$t /1000
allAverageDat_ica_exp1_ps_n %>% filter(t <= 0.65 & t >= 0.25) %>%
  group_by(SubName) %>% nest() %>%
  mutate(model = map(data, cnv_fit)) %>% 
  unnest(model, .drop = TRUE) -> lm_exp1_ps_slope1
lm_exp1_ps_slope1$cond <- "PS"

t.test(lm_exp1_ps_slope1$slope, mu = 0, alternative = "two.sided")
```

```{r}
#based on condition?
lm_exp1_slope1 <- rbind(lm_exp1_ps_slope1, lm_exp1_ns_slope1)

lm_exp1_slope1 %>% group_by(cond) %>%dplyr::summarise(n=n(), m_slope=mean(slope), se_slope=sd(slope)/sqrt(n-1))
```

```{r}
ezANOVA(data = lm_exp1_slope1, dv= slope, wid=SubName, within=.(cond))
```

# 1.1 CNV activity

```{r}
m_allAverageDat_exp1 <- allAverageDat_exp1 %>% group_by(targetDur,t, cond) %>% 
  dplyr::summarise(n= n(), m_CNV = mean(CNV),se_CNV = sd(CNV)/sqrt(n-1))
```


###separate plots
```{r}
c1 <- c("#000333", "#0000CC", "#3399CC",  "#33FFFF","#66CC00", "#CC0066", "#FF3333")
p1 <- ggplot(m_allAverageDat_exp1 %>% filter(cond=="PS"), aes(x=t, y=m_CNV, color=as.factor(targetDur)))+geom_line(size=0.5)+
  scale_color_manual(values=c1) + mytheme+
  scale_y_reverse() + labs(x = "Time (ms)", y = "Amplitude (\u03BCV)", color='') +
  scale_x_continuous(breaks= c(0, 400,800,1200, 1600, 2000), labels = c(0, 400,800,1200, 1600,2000)) +coord_cartesian(ylim = c(7.5, -6),xlim = c(-200, 2400), expand = c(0,0))+ ggtitle("Short context (PS)")+geom_hline(yintercept=0, linetype="dashed")+theme(legend.key.size = unit(0.5, "cm"))

p1<-p1+theme(axis.title.y = element_text(angle=90)) 
p1
```

```{r}
c11<-c("#000333", "#3833cc", "#33a1cc","#33cca6","#c2cc33", "#b233cc", "#FF3333")

p11<-ggplot(m_allAverageDat_exp1 %>% filter(cond=="NS" ), aes(x=t, y=m_CNV, color=as.factor(targetDur)))+geom_line(size=0.5)+
  scale_color_manual(values=c11)+
  scale_y_reverse() + mytheme+ labs(x = "Time (ms)", y = "Amplitude (\u03BCV)", color='') +
  scale_x_continuous(breaks= c(0, 400,800,1200, 1600,2000), labels = c(0, 400,800,1200, 1600,2000)) +coord_cartesian(ylim = c(7.5, -6),xlim = c(-200, 2400), expand = c(0,0))+ ggtitle("Long context (NS)")+theme(legend.key.size = unit(0.5, "cm"))+
  geom_hline(yintercept=0, linetype="dashed")

p11<-p11 +theme(axis.title.y = element_text(angle=90))
p11
```


# 1.2 CNV-Peak

We check the negativity peak latency within 0-1600 ms using a custom peak detection function.

```{r}
mexp1_negPeak_lat <- exp1_negPeak_lat %>% group_by(cond, SubName) %>% dplyr::summarise(m_lat = mean(Latency), m_amp= mean(Amplitude))

mmexp1_negPeak_lat <-mexp1_negPeak_lat %>% group_by(cond) %>% dplyr::summarise(n=n(),mm_lat = mean(m_lat), se_lat= sd(m_lat)/sqrt(n-1), mm_amp=mean(m_amp), se_amp=sd(m_amp)/sqrt(n-1))
mmexp1_negPeak_lat
```

```{r}
ezANOVA(data = mexp1_negPeak_lat, dv= m_lat, wid=SubName, within=.(cond))
```

```{r}
ezANOVA(data = mexp1_negPeak_lat, dv= m_amp, wid=SubName, within=.(cond))
```



##barplot
```{r}
exp1_lat= ggplot(data=mmexp1_negPeak_lat%>% mutate(cond = factor(cond, levels = c("PS", "NS"))), aes(x=cond, y=mm_lat, color=cond,fill=cond)) +geom_bar(stat="identity",position=position_dodge(), width=0.6, fill="white", alpha=0.9)+  mytheme+ scale_color_manual(values=colors_plot)+ theme(legend.position = "") + labs(x="", y="Peak latency (ms)") + geom_errorbar(data=mmexp1_negPeak_lat%>% mutate(cond = factor(cond, levels = c("PS", "NS"))),aes(ymin = mm_lat-se_lat,ymax = mm_lat+ se_lat),width = 0.1, position = position_dodge(width=0.6))+ coord_cartesian(ylim=c(0,1000)) +theme(axis.text.y= element_text(angle=90))+ggtitle("")+theme(plot.title = element_text(color = "black", size = 10))
exp1_lat
```

```{r}
exp1_amp= ggplot(data=mmexp1_negPeak_lat%>% mutate(cond = factor(cond, levels = c("PS", "NS"))), aes(x=cond, y=mm_amp, color=cond,fill=cond)) +geom_bar(stat="identity",position=position_dodge(), width=0.6, fill="white", alpha=0.9)+ scale_y_reverse() +mytheme+ scale_color_manual(values=colors_plot)+ theme(legend.position = "") + labs(x="", y="Peak amplitude (\u03BCV)") + geom_errorbar(aes(ymin = mm_amp-se_amp,ymax = mm_amp+ se_amp),width = 0.1, position = position_dodge(width=0.6))+coord_cartesian(ylim=c(0,-7))+ggtitle("")+theme(plot.title = element_text(color = "black", size = 10))
exp1_amp
```

```{r}
m_peak_exp1_targets <- exp1_negPeak_lat %>% group_by(targetDur,cond) %>% dplyr::summarise(n=n(), m_lat = mean(Latency),
                   se_lat = sd(Latency)/sqrt(n-1), 
                   m_amp = mean(Amplitude), 
                   se_amp = sd(Amplitude)/sqrt(n-1))
```


##plot
```{r}
p_exp1_lat <-  m_peak_exp1_targets%>% mutate(cond = factor(cond, levels = c("PS", "NS")))%>%
  ggplot(aes(x=targetDur, y=m_lat, color=cond))+geom_point()+ 
  geom_abline(slope=1, linetype="dashed", color='black')+
  geom_errorbar(aes(ymin = m_lat-se_lat,ymax = m_lat+se_lat),width = 40, position = position_dodge(width=0.5))+mytheme+
   theme(axis.text.y= element_text(angle=90)) +
  geom_line()+theme(legend.position = c(0.9,0.9))+
    labs(x = "Target interval (ms)", y = "Peak latency (ms)") +scale_color_manual(values = colors_plot) + theme(legend.title = element_blank())+ #scale_x_continuous(breaks=c(400,800,1200,1600), labels = c("0.4","0.8","1.2","1.6"))
  coord_cartesian(ylim = c(500,1200)) +theme(axis.text.y= element_text(angle=90))+ggtitle("CNV")

 p_exp1_lat
```

#1.3 CNV-mean amplitude
```{r}
# calculate average CNV
t_ps <- data.frame()
targetList <- c(400,504,636,800,1008,1270, 1600)
for (target in targetList){
  m_cnv_1 <- allAverageDat_exp1%>% filter( (t>= 250 & t <= (target+250))) %>% group_by(cond,targetDur,SubName) %>% dplyr::summarise(n= n(), m_CNV = mean(CNV),se_CNV = sd(CNV)/sqrt(n-1))
  m_cnv_1= m_cnv_1%>% filter(targetDur== target & cond=="PS")
  t_ps <- rbind2(t_ps,m_cnv_1, by=c("SubName","targetDur","cond")) 
}
```

```{r}
# calculate average CNV
t_ns <- data.frame()
targetList <- c(400,730,992,1200,1366,1496, 1600)
for (target in targetList){
  m_cnv_2 <- allAverageDat_exp1%>% filter( (t>= 250 & t <= (target+250))) %>% group_by(cond,targetDur,SubName) %>% dplyr::summarise(n= n(), m_CNV = mean(CNV),se_CNV = sd(CNV)/sqrt(n-1))
  m_cnv_2= m_cnv_2%>% filter(targetDur== target & cond=="NS")
  t_ns <- rbind2(t_ns,m_cnv_2, by=c("SubName","targetDur","cond")) 
}
```



```{r}
exp1_t <- rbind(t_ps, t_ns)
t_targets <- exp1_t %>% group_by(targetDur, cond)%>% dplyr::summarise(n=n(), mm_CNV =mean(m_CNV), se_CNV=sd(m_CNV)/sqrt(n-1))

m_exp1_t <- exp1_t %>% group_by(SubName, cond)%>% dplyr::summarise(n=n(), mm_CNV =mean(m_CNV), se_CNV=sd(m_CNV)/sqrt(n-1))

m_exp1_t %>% group_by(cond)%>% dplyr::summarise(n=n(), mmm_CNV =mean(mm_CNV))
```

##plot

```{r}
exp1_ave <-  t_targets%>% mutate(cond = factor(cond, levels = c("PS", "NS")))%>% 
  ggplot(aes(x=targetDur, y=mm_CNV, color=cond)) + 
  geom_point(position = position_dodge(width=30))+
  geom_errorbar(aes(ymin = mm_CNV-se_CNV,ymax = mm_CNV+se_CNV),width = 80, position = position_dodge(width=30),size=0.5)+mytheme+
  geom_line()+
    labs(x = "Target interval (ms)", y = "Mean amplitude (\u03BCV)") +scale_y_reverse()+
    scale_color_manual(values = colors_plot)+ theme(legend.position = c(0.9,0.25), legend.title = element_blank()) + scale_x_continuous(breaks = c(400,800,1200,1600), labels = c(400,800,1200,1600))+coord_cartesian(ylim = c(2.5, -3)) + ggtitle("")

 exp1_ave
```
##mixed linear model
```{r}
exp1_t$targetDur = exp1_t$targetDur /1000
exp1_t$cond = as.factor(exp1_t$cond)
contrasts(exp1_t$cond) = contr.Sum(levels(exp1_t$cond))
```

```{r}
# fit the linear mixed model
exp1_mean_cnv = lmer(m_CNV ~
              cond*targetDur + # fixed effect, covariate and interaction
              (1+targetDur|SubName), # random effect
              data = exp1_t)
#output as table
tab_model(exp1_mean_cnv, p.val = 'kr') # Kenward-Roger approximation, alternative 'satterthwaite' approxmiation
```

###Figure 3

```{r}
fig2_11 <- cowplot::plot_grid(p1,exp1_lat, exp1_amp,
                               nrow = 1, labels = c("a", "b","c"), rel_widths = c(2.5,0.5,0.5))
fig2_12 <- cowplot::plot_grid(p11,exp1_ave,
                               nrow = 1, labels = c("d", "e"), rel_widths = c(2.5,1))
fig2_13 <- cowplot::plot_grid(fig2_11, fig2_12, nrow=2, rel_heights = c(1,1))

ggsave("figures/figure3.png", fig2_13, width = 8, height = 4.5) 
fig2_13       
```




# 2.Offset-P2 

```{r}
m_data_ica_p2_exp1 <- data_ica_p2_exp1 %>% group_by(targetDur,t, cond) %>% 
  dplyr::summarise(n= n(), m_P2 = mean(P2),
                   se_P2 = sd(P2)/sqrt(n-1))

m_data_ica_p2_exp1$targetDur <- as.factor(m_data_ica_p2_exp1$targetDur)
```

```{r}
data_ica_p2_exp1_sub= data_ica_p2_exp1 %>% group_by(SubName, targetDur,t, cond) %>% 
  dplyr::summarise(n= n(), m_P2 = mean(P2),
                   se_P2 = sd(P2)/sqrt(n-1))
```

##separate plots
```{r}
colors_4_exp1 <- c("#000333", "#3833cc", "#33a1cc","#33cca6","#c2cc33", "#b233cc", "#FF3333")

plot_4_p2_exp1_ps <- ggplot(m_data_ica_p2_exp1 %>% filter(cond=="PS"), aes(x=t, y=m_P2, color=as.factor(targetDur)))+geom_line(size=0.5)+
  scale_color_manual(values=colors_4_exp1)+mytheme+
  scale_y_reverse() +labs(x = "Time (ms)", y = "Amplitude (\u03BCV)", color='') +
  scale_x_continuous(breaks= c(0, 200,400,600), labels = c(0, 200,400,600)) +coord_cartesian(ylim = c(10,-5),xlim = c(-100, 700), expand = c(0,0))+ ggtitle("Short context (PS)")+theme(legend.position = "right")+theme(axis.title.y = element_text(angle=90))+geom_hline(yintercept=0, linetype="dashed")+theme(legend.key.size = unit(0.5, "cm"))+
  annotate("rect", xmin = 140, xmax = 300, ymin = -3, ymax = 9,
           alpha = 0.1,fill = "gray", color="darkgray")+annotate("text", x = 225, y = -4, label = "P2")

plot_4_p2_exp1_ps
```


```{r}
colors_4_exp1 <- c("#000333", "#0000CC", "#3399CC","#33FFFF","#66CC00", "#CC0066", "#FF3333")

plot_4_p2_exp1_ns <- ggplot(m_data_ica_p2_exp1 %>% filter(cond=="NS"), aes(x=t, y=m_P2, color=as.factor(targetDur)))+geom_line(size=0.5)+
  scale_color_manual(values=colors_4_exp1)+mytheme+
  scale_y_reverse() + labs(x = "Time (ms)", y = "Amplitude (\u03BCV)", color='') +
  scale_x_continuous(breaks= c(0, 200,400,600), labels = c(0, 200,400,600)) +coord_cartesian(ylim = c(10,-5),xlim = c(-100, 700), expand = c(0,0))+ ggtitle("Long context (NS)")+theme(legend.position = "right")+theme(axis.title.y = element_text(angle=90))+geom_hline(yintercept=0, linetype="dashed")+theme(legend.key.size = unit(0.5, "cm"))+
  annotate("rect", xmin = 140, xmax = 300, ymin = -3, ymax = 9,
           alpha = 0.1,fill = "gray", color="darkgray")+annotate("text", x = 225, y = -4, label = "P2")


plot_4_p2_exp1_ns
```



#2.1Offset-P2-Peak
```{r}
m_exp1_p2_peaks <- exp1_pos_peaks %>% group_by(cond, SubName) %>%
  dplyr::summarise(n=n(), m_lat =mean(Latency), m_amp = mean(Amplitude))

m_exp1_p2_peaks_targets <- exp1_pos_peaks %>% group_by(cond, targetDur) %>%
  dplyr::summarise(n=n(), m_lat =mean(Latency), m_amp = mean(Amplitude),se_lat= sd(Latency)/ sqrt(n-1), se_amp = sd(Amplitude)/ sqrt(n-1))

mm_exp1_p2_peaks <- m_exp1_p2_peaks %>% group_by(cond) %>%
  dplyr::summarise(n=n(), mm_lat =mean(m_lat), mm_amp = mean(m_amp), se_lat= sd(m_lat)/ sqrt(n-1), se_amp = sd(m_amp)/ sqrt(n-1))
mm_exp1_p2_peaks
```
###mixed linear model
```{r}
exp1_pos_peaks$targetDur = exp1_pos_peaks$targetDur /1000
#exp1_pos_peaks$Latency = exp1_pos_peaks$Latency /1000

exp1_pos_peaks$cond = as.factor(exp1_pos_peaks$cond)
contrasts(exp1_pos_peaks$cond) = contr.Sum(levels(exp1_pos_peaks$cond))
```

```{r}
# fit the linear mixed model
mexp1_p2_peaks = lmer(Latency ~
              cond*targetDur + # fixed effect, covariate and interaction
              (1+targetDur|SubName), # random effect
              data = exp1_pos_peaks)
#output as table
tab_model(mexp1_p2_peaks, p.val = 'kr') # Kenward-Roger approximation, alternative 'satterthwaite' approxmiation
```


##plot

```{r}
p2p3_exp1_lat <-  m_exp1_p2_peaks_targets%>% mutate(cond = factor(cond, levels = c("PS", "NS")))%>% 
  ggplot(aes(x=targetDur, y=m_lat, color=cond))+geom_point()+
  geom_errorbar(aes(ymin = m_lat-se_lat,ymax = m_lat+se_lat),width = 40, position = position_dodge(width=0.5))+mytheme+
    labs(x = "Target interval (ms)", y = "Peak latency (ms)") +
    scale_color_manual(values = colors_plot)+ theme(legend.position = c(0.9,0.9), legend.title = element_blank())+ scale_x_continuous(breaks=c(400,800,1200,1600)) + geom_line() +ggtitle("P2")
#+ coord_cartesian(ylim = c(-3, -8)) 

 p2p3_exp1_lat
```

#2.2 P2-mean amplitude

```{r}
p2_amp_exp1 <- data_ica_p2_exp1 %>%filter(between(t, 140, 300))%>% group_by(SubName, cond, targetDur) %>% dplyr::summarise(n=n(), m_p2=mean(P2), se_p2= sd(P2)/sqrt(n-1))
```

## mixed linear model
```{r}
p2_amp_exp1$targetDur = p2_amp_exp1$targetDur /1000
p2_amp_exp1$cond = as.factor(p2_amp_exp1$cond)
contrasts(p2_amp_exp1$cond) = contr.Sum(levels(p2_amp_exp1$cond))
```


```{r}
# fit the linear mixed model
mlm = lmer(m_p2 ~
              cond*targetDur + # fixed effect, covariate and interaction
              (1+targetDur|SubName), # random effect
              data = p2_amp_exp1)
#output as table
tab_model(mlm, p.val = 'kr') # Kenward-Roger approximation, alternative 'satterthwaite' approxmiation
```

##plot
```{r}
p2_amp_exp1 <- data_ica_p2_exp1 %>%filter(between(t, 140, 300))%>% group_by(SubName, cond, targetDur) %>% dplyr::summarise(n=n(), m_p2=mean(P2), se_p2= sd(P2)/sqrt(n-1))
p2_amp_exp1_t <- p2_amp_exp1 %>%  group_by(cond, targetDur) %>% dplyr::summarise(n=n(), mm_p2=mean(m_p2), se_p2= sd(m_p2)/sqrt(n-1))
p2_amp_exp1_tt <- p2_amp_exp1_t %>%  group_by(cond) %>% dplyr::summarise(n=n(), mmm_p2=mean(mm_p2))
p2_amp_exp1_tt
```

```{r}
p2p3_exp1_amp_mean <-  p2_amp_exp1_t%>% mutate(cond = factor(cond, levels = c("PS", "NS")))%>% 
  ggplot(aes(x=targetDur, y=mm_p2, color=cond))+geom_point()+
  geom_errorbar(aes(ymin = mm_p2-se_p2,ymax = mm_p2+se_p2),width = 40, position = position_dodge(width=0.5))+mytheme+
    labs(x = "Target interval (ms)", y = "Mean amplitude (\u03BCV)") +
    scale_color_manual(values = colors_plot)+ theme(legend.position = "", legend.title = element_blank())+ scale_x_continuous(breaks=c(400,800,1200,1600)) + geom_line() + coord_cartesian(ylim = c(-2, 9))+ theme(legend.position = c(0.9,0.24), legend.title = element_blank()) +ggtitle("P2")

 p2p3_exp1_amp_mean
```


###Figure 4

```{r}
fig_exp1_p2_ns_ps <- cowplot::plot_grid(plot_4_p2_exp1_ps,p2p3_exp1_lat, plot_4_p2_exp1_ns, p2p3_exp1_amp_mean,
                               nrow = 2, labels = c("a", "b","c","d"), rel_widths = c(3,1.3))

ggsave("figures/figure4.png", fig_exp1_p2_ns_ps, width = 8, height = 4.5) 

fig_exp1_p2_ns_ps
```


#EXPERIMENT 2

```{r}
allAverageDat_ica_exp2_af <- allAverageDat_exp2 %>% filter(cond== "AF")
allAverageDat_ica_exp2_df <- allAverageDat_exp2 %>% filter(cond== "DF")
```

# 1.CNV

```{r}
m_allAverageDat_exp2 <- allAverageDat_exp2 %>% group_by(targetDur,t, cond) %>% 
  dplyr::summarise(n= n(), m_CNV = mean(CNV),se_CNV = sd(CNV)/sqrt(n-1))
```

## separate plots 
```{r}
c2 <- c("#000333", "#0000CC","#3399CC", "#33FFFF","#FF9933", "#CC0066", "#FF3333")
p2 <-ggplot(m_allAverageDat_exp2 %>% filter(cond=="DF" ), aes(x=t, y=m_CNV, color=as.factor(targetDur)))+geom_line(size=0.5 )+
  scale_color_manual(values=c2)+ mytheme+
  scale_y_reverse() +
  theme(axis.ticks.length=unit(.10, "cm"))+ labs(x = "Time (ms)", y = "Amplitude (\u03BCV)", color='') +
  scale_x_continuous(breaks= c(0, 400,800,1200, 1600,2000), labels = c(0, 400,800,1200, 1600,2000)) +coord_cartesian(ylim = c(7.5, -6),xlim = c(-200, 2400), expand = c(0,0))+ ggtitle("Short context (DF)")+ theme(legend.position = "right")+ geom_hline(yintercept = 0, linetype="dashed")+theme(legend.key.size = unit(0.5, "cm"))

p2<-p2 +theme(axis.title.y = element_text(angle=90))
p2
```



```{r}
c22 <-c("#000333", "#0000CC", "#3399CC","#33FFFF","#FF9933", "#CC0066", "#FF3333")
p22<-ggplot(m_allAverageDat_exp2 %>% filter(cond=="AF" ), aes(x=t, y=m_CNV, color=as.factor(targetDur)))+geom_line(size=0.5)+
  scale_color_manual(values=c22)+mytheme+
  scale_y_reverse() +
  theme(axis.ticks.length=unit(.10, "cm"))+ labs(x = "Time (ms)", y = "Amplitude (\u03BCV)", color='') +
  scale_x_continuous(breaks= c(0, 400,800,1200, 1600,2000), labels = c(0, 400,800,1200, 1600,2000)) +coord_cartesian(ylim = c(7.5, -6),xlim = c(-200, 2400), expand = c(0,0))+ ggtitle("Long context (AF)")+theme(legend.position = "right")+ geom_hline(yintercept = 0, linetype="dashed")+theme(legend.key.size = unit(0.5, "cm"))

p22<-p22 +theme(axis.title.y = element_text(angle=90))
p22
```


## Slopes

```{r}
allAverageDat_ica_exp2_af_n = allAverageDat_ica_exp2_af
allAverageDat_ica_exp2_af_n$targetDur = allAverageDat_ica_exp2_af_n$targetDur /1000
allAverageDat_ica_exp2_af_n$t = allAverageDat_ica_exp2_af_n$t /1000
allAverageDat_ica_exp2_af_n %>% filter(t <= 0.65 & t >= 0.25) %>% 
  group_by(SubName) %>% nest() %>%
  mutate(model = map(data, cnv_fit)) %>% 
  unnest(model, .drop = TRUE) -> lm_exp2_af_slope1
lm_exp2_af_slope1$cond <- "AF"

t.test(lm_exp2_af_slope1$slope, mu = 0, alternative = "two.sided")
```


```{r}
#in DF
allAverageDat_ica_exp2_df_n = allAverageDat_ica_exp2_df
allAverageDat_ica_exp2_df_n$targetDur = allAverageDat_ica_exp2_df_n$targetDur /1000
allAverageDat_ica_exp2_df_n$t = allAverageDat_ica_exp2_df_n$t /1000
allAverageDat_ica_exp2_df_n %>% filter(t <= 0.65 & t >= 0.25) %>%
  group_by(SubName) %>% nest() %>%
  mutate(model = map(data, cnv_fit)) %>% 
  unnest(model, .drop = TRUE) -> lm_exp2_df_slope1
lm_exp2_df_slope1$cond <- "DF"

t.test(lm_exp2_df_slope1$slope, mu = 0, alternative = "two.sided")
```

```{r}
#based on condition?
lm_exp2_slope1 <- rbind(lm_exp2_df_slope1, lm_exp2_af_slope1)

lm_exp2_slope1 %>% group_by(cond) %>%dplyr::summarise(n=n(), m_slope=mean(slope), se_slope=sd(slope)/sqrt(n-1))
```

```{r}
#Do they differ in different conditions?
ezANOVA(data = lm_exp2_slope1, dv= slope, wid=SubName, within=.(cond))
```



# 1.2.CNV-Peak
```{r}
mexp2_negPeak_lat <- exp2_negPeak_lat %>% group_by(cond, SubName) %>% dplyr::summarise(m_lat = mean(Latency), m_amp= mean(Amplitude))

mmexp2_negPeak_lat <- mexp2_negPeak_lat %>% group_by(cond) %>% dplyr::summarise(n=n(),mm_lat = mean(m_lat), se_lat= sd(m_lat)/sqrt(n-1), mm_amp=mean(m_amp), se_amp= sd(m_amp)/sqrt(n-1))
mmexp2_negPeak_lat
```
```{r}
ezANOVA(data = mexp2_negPeak_lat, dv= m_lat, wid=SubName, within=.(cond))
```

```{r}
ezANOVA(data = mexp2_negPeak_lat, dv= m_amp, wid=SubName, within=.(cond))
```



```{r}
m_peak_exp2_targets <- exp2_negPeak_lat %>% group_by(targetDur,cond) %>% dplyr::summarise(n=n(), m_lat = mean(Latency),
                   se_lat = sd(Latency)/sqrt(n-1), 
                   m_amp = mean(Amplitude), 
                   se_amp = sd(Amplitude)/sqrt(n-1))

m_peak_exp2_subs <- exp2_negPeak_lat %>% group_by(SubName,cond) %>% dplyr::summarise(n=n(), m_lat = mean(Latency),
                   se_lat = sd(Latency)/sqrt(n-1), 
                   m_amp = mean(Amplitude), 
                   se_amp = sd(Amplitude)/sqrt(n-1))
```

##barplot
```{r}
exp2_lat= ggplot(data=mmexp2_negPeak_lat%>% mutate(cond = factor(cond, levels = c("DF", "AF"))), aes(x=cond, y=mm_lat, color=cond,fill=cond)) +geom_bar(stat="identity",position=position_dodge(), width=0.6, fill="white", alpha=0.9)+  mytheme+ scale_color_manual(values=colors_plot)+ theme(legend.position = "") + labs(x="", y="Peak latency (ms)") + geom_errorbar(aes(ymin = mm_lat-se_lat,ymax = mm_lat+ se_lat),width = 0.1, position = position_dodge(width=0.6)) +coord_cartesian(ylim=c(0,1000))+ ggtitle("")+theme(axis.text.y= element_text(angle=90))+theme(plot.title = element_text(color = "black", size = 10))
exp2_lat
```

```{r}
exp2_amp= ggplot(data=mmexp2_negPeak_lat%>% mutate(cond = factor(cond, levels = c("DF", "AF"))), aes(x=cond, y=mm_amp, color=cond,fill=cond)) +geom_bar(stat="identity",position=position_dodge(), width=0.6, fill="white", alpha=0.9)+ scale_y_reverse() +mytheme+ scale_color_manual(values=colors_plot)+ theme(legend.position = "") + labs(x="", y="Peak amplitude (\u03BCV)") + geom_errorbar(aes(ymin = mm_amp-se_amp,ymax = mm_amp+ se_amp),width = 0.1, position = position_dodge(width=0.6))+coord_cartesian(ylim=c(0,-7))+ggtitle("")+theme(plot.title = element_text(color = "black", size = 10))
exp2_amp
```

##plot
```{r}
p_exp2_lat <-  m_peak_exp2_targets%>% mutate(cond = factor(cond, levels = c("DF", "AF")))%>% 
  ggplot(aes(x=targetDur, y=m_lat, color=cond))+ geom_point(position = position_dodge(width=30))+
  geom_errorbar(aes(ymin = m_lat-se_lat,ymax = m_lat+se_lat),width = 30, position = position_dodge(width=30))+mytheme+ geom_abline(slope=1, linetype="dashed", color='black')+
    labs(x = "Target interval (ms)", y = "Peak latency (ms)") +
   scale_color_manual(values = colors_plot)+ theme(legend.position = c(0.9,0.24), legend.title = element_blank())+ 
  geom_line()+coord_cartesian(ylim = c(500,1200))+theme(axis.text.y= element_text(angle=90)) +ggtitle("CNV")

 p_exp2_lat
```


# 1.3 CNV-mean amplitude
```{r}
# calculate average CNV
t_df <- data.frame()
targetList <- c(400,600,800,1000,1200,1400, 1600)
for (target in targetList){
  m_cnv_1 <- allAverageDat_exp2%>% filter( (t>= 250 & t <= (target+250))) %>% group_by(cond,targetDur,SubName) %>% dplyr::summarise(n= n(), m_CNV = mean(CNV),se_CNV = sd(CNV)/sqrt(n-1))
  m_cnv_1= m_cnv_1%>% filter(targetDur== target & cond=="DF")
  t_df <- rbind2(t_df,m_cnv_1, by=c("SubName","targetDur","cond")) 
}
```


```{r}
# calculate average CNV
t_af <- data.frame()
targetList <- c(400,600,800,1000,1200,1400, 1600)
for (target in targetList){
  m_cnv_2 <- allAverageDat_exp2%>% filter( (t>= 250 & t <= (target+250))) %>% group_by(cond,targetDur,SubName) %>% dplyr::summarise(n= n(), m_CNV = mean(CNV),se_CNV = sd(CNV)/sqrt(n-1))
  m_cnv_2= m_cnv_2%>% filter(targetDur== target & cond=="AF")
  t_af <- rbind2(t_af,m_cnv_2, by=c("SubName","targetDur","cond")) 
}
```


```{r}
t_exp2 <- rbind(t_df, t_af)
t_targets_exp2 <- t_exp2 %>% group_by(targetDur, cond)%>% dplyr::summarise(n=n(), mm_CNV =mean(m_CNV), se_CNV=sd(m_CNV)/sqrt(n-1))

m_t_exp2 <- t_exp2 %>% group_by(SubName, cond)%>% dplyr::summarise(n=n(), mm_CNV =mean(m_CNV), se_CNV=sd(m_CNV)/sqrt(n-1))
m_t_exp2 %>% group_by(cond)%>% dplyr::summarise(n=n(), mmm_CNV =mean(mm_CNV))
```


##plot

```{r}
exp2_ave <-  t_targets_exp2%>% mutate(cond = factor(cond, levels = c("DF", "AF")))%>% 
  ggplot(aes(x=targetDur, y=mm_CNV, color=cond)) + 
  geom_point(position = position_dodge(width=30))+
  geom_errorbar(aes(ymin = mm_CNV-se_CNV,ymax = mm_CNV+se_CNV),width = 80, position = position_dodge(width=30),size=0.5)+mytheme+
  geom_line()+
    labs(x = "Target intervals (ms)", y = "Mean amplitude (\u03BCV)") +scale_y_reverse()+
    scale_color_manual(values = colors_plot)+ theme(legend.position = c(0.9,0.24), legend.title = element_blank()) + scale_x_continuous(breaks = c(400,800,1200,1600), labels = c(400,800,1200,1600))+coord_cartesian(ylim = c(2.5, -3))+ ggtitle("")

 exp2_ave
```


###mixed linear model

```{r}
t_exp2$targetDur = t_exp2$targetDur /1000
```

```{r}
t_exp2$cond = as.factor(t_exp2$cond)
contrasts(t_exp2$cond) = contr.Sum(levels(t_exp2$cond))

# fit the linear mixed model
m_exp1_cnv_latency = lmer(m_CNV ~
              cond*targetDur + # fixed effect, covariate and interaction
              (1+targetDur|SubName), # random effect
              data = t_exp2)
#output as table
tab_model(m_exp1_cnv_latency, p.val = 'kr') # Kenward-Roger approximation, alternative 'satterthwaite' approxmiation
```


###Figure 7

```{r}
fig5_11 <- cowplot::plot_grid(p2,exp2_lat, exp2_amp,
                               nrow = 1, labels = c("a", "b","c"), rel_widths = c(2.5,0.5,0.5))
fig5_12 <- cowplot::plot_grid(p22,exp2_ave,
                               nrow = 1, labels = c("d", "e"), rel_widths = c(2.5,1))
fig5_13 <- cowplot::plot_grid(fig5_11, fig5_12, nrow=2, rel_heights = c(1,1))

ggsave("figures/figure7.png", fig5_13, width = 8, height = 4.5) 
fig5_13       
```


# 2. LPCt
```{r}
m_data_ica_p2 <- data_ica_pos_exp2 %>% group_by(targetDur,t, cond) %>% 
  dplyr::summarise(n= n(), m_P2 = mean(P2),
                   se_P2 = sd(P2)/sqrt(n-1))
```


## separate plots
```{r}
colors_1_exp2 <- c("#000333", "#0000CC","#3399CC", "#33FFFF","#FF9933", "#CC0066", "#FF3333")

plot_1_p2_exp2_df <- ggplot(m_data_ica_p2 %>% filter(cond=="DF"), aes(x=t, y=m_P2, color=as.factor(targetDur)))+geom_line(size=0.5)+
  scale_color_manual(values=colors_1_exp2)+ mytheme+
  scale_y_reverse() +labs(x = "Time (ms)", y = "Amplitude (\u03BCV)", color='') +
  #scale_x_continuous(breaks= c(0, 200,400,600), labels = c(0, 200,400,600)) +
  coord_cartesian(ylim = c(10,-5), expand = c(0,0))+ ggtitle("Short context (DF)")+theme(legend.position = "right")+theme(axis.title.y = element_text(angle=90)) +theme(legend.position = "right")+theme(axis.title.y = element_text(angle=90))+geom_hline(yintercept=0, linetype="dashed")+theme(legend.key.size = unit(0.5, "cm"))+
  annotate("rect", xmin = 300, xmax = 500, ymin = -3, ymax = 9,
           alpha = 0.1,fill = "gray", color="darkgray")+annotate("text", x = 400, y = -4, label = "LPCt")

plot_1_p2_exp2_df
```


```{r}
colors_1_exp2 <- c("#000333", "#0000CC","#3399CC", "#33FFFF","#FF9933", "#CC0066", "#FF3333")

plot_1_p2_exp2_af <- ggplot(m_data_ica_p2 %>% filter(cond=="AF"), aes(x=t, y=m_P2, color=as.factor(targetDur)))+geom_line(size=0.5)+
  scale_color_manual(values=colors_1_exp2)+ mytheme+
  scale_y_reverse() +labs(x = "Time (ms)", y = "Amplitude (\u03BCV)", color='') +
  scale_x_continuous(breaks= c(0, 200,400,600), labels = c(0, 200,400,600)) +coord_cartesian(ylim = c(10,-5),xlim = c(-100, 700), expand = c(0,0))+ ggtitle("Long context (AF)")+theme(legend.position = "right")+theme(legend.key.size = unit(0.5, "cm"))+theme(axis.title.y = element_text(angle=90))+ geom_hline(yintercept = 0, linetype="dashed")+ 
  annotate("rect", xmin = 300, xmax = 500, ymin = -3, ymax = 9,
           alpha = 0.1,fill = "gray", color="darkgray")+annotate("text", x = 400, y = -4, label = "LPCt")

plot_1_p2_exp2_af
```



# 2.1 LPCt-Peak
```{r}
m_exp2_p2_peaks <- exp2_pos_peaks %>% group_by(cond, SubName) %>%
  dplyr::summarise(n=n(), m_lat =mean(Latency), m_amp = mean(Amplitude))

m_exp2_p2_peaks_targets <- exp2_pos_peaks %>% group_by(cond, targetDur) %>%
  dplyr::summarise(n=n(), m_lat =mean(Latency), m_amp = mean(Amplitude), se_lat= sd(Latency)/ sqrt(n-1), se_amp = sd(Amplitude)/ sqrt(n-1))

mm_exp2_p2_peaks <- m_exp2_p2_peaks %>% group_by(cond) %>%
  dplyr::summarise(n=n(), mm_lat =mean(m_lat), mm_amp = mean(m_amp), se_lat= sd(m_lat)/ sqrt(n-1), se_amp = sd(m_amp)/ sqrt(n-1))
mm_exp2_p2_peaks
```

###mixed linear model
```{r}
exp2_pos_peaks$targetDur = exp2_pos_peaks$targetDur /1000
#exp2_pos_peaks$Latency = exp2_pos_peaks$Latency /1000
exp2_pos_peaks$cond = as.factor(exp2_pos_peaks$cond)
contrasts(exp2_pos_peaks$cond) = contr.Sum(levels(exp2_pos_peaks$cond))
```

```{r}
# fit the linear mixed model
mexp2_p2_peaks = lmer(Latency ~
              cond*targetDur + # fixed effect, covariate and interaction
              (1+targetDur|SubName), # random effect
              data = exp2_pos_peaks)
#output as table
tab_model(mexp2_p2_peaks, p.val = 'kr') # Kenward-Roger approximation, alternative 'satterthwaite' approxmiation
```


##plot
```{r}
p2p3_exp2_lat <-  m_exp2_p2_peaks_targets%>% mutate(cond = factor(cond, levels = c("DF", "AF")))%>% 
  ggplot(aes(x=targetDur, y=m_lat, color=cond))+geom_point()+
  geom_errorbar(aes(ymin = m_lat-se_lat,ymax = m_lat+se_lat),width = 40, position = position_dodge(width=0.5))+mytheme+
    labs(x = "Target interval (ms)", y = "Peak latency (ms)") +
    scale_color_manual(values = colors_plot)+ theme(legend.position = c(0.9,0.9), legend.title = element_blank())+ scale_x_continuous(breaks=c(400,800,1200,1600)) + geom_line() +ggtitle("LPCt")
#+ coord_cartesian(ylim = c(-3, -8))

 p2p3_exp2_lat
```

#2.2- mean amplitude

```{r}
# test target Dur x condition x mean amplitude
mdata_ica_p2_exp2 <- data_ica_pos_exp2%>% filter( (t>=300 & t<=500)) %>% group_by(targetDur,cond, SubName) %>% dplyr::summarise(n= n(), m_pos = mean(P2),se_pos = sd(P2)/sqrt(n-1))
```


###mixed linear model
```{r}
mdata_ica_p2_exp2$targetDur = mdata_ica_p2_exp2$targetDur /1000
```

```{r}
mdata_ica_p2_exp2$cond = as.factor(mdata_ica_p2_exp2$cond)
contrasts(mdata_ica_p2_exp2$cond) = contr.Sum(levels(mdata_ica_p2_exp2$cond))

# fit the linear mixed model
m_lm = lmer(m_pos ~
              cond*targetDur + # fixed effect, covariate and interaction
              (1+targetDur|SubName), # random effect
              data = mdata_ica_p2_exp2)
#output as table
tab_model(m_lm, p.val = 'kr') # Kenward-Roger approximation, alternative 'satterthwaite' approxmiation
```


##plot
```{r}
mdata_ica_p2_exp2 <- data_ica_pos_exp2%>% filter( (t>=300 & t<=500)) %>% group_by(targetDur,cond, SubName) %>% dplyr::summarise(n= n(), m_P2 = mean(P2),se_P2 = sd(P2)/sqrt(n-1))

mmdata_ica_p2_exp2 <- mdata_ica_p2_exp2 %>% group_by(targetDur,cond) %>% dplyr::summarise(n= n(), mm_P2 = mean(m_P2),se_P2 = sd(m_P2)/sqrt(n-1))

m_p2_exp2 <-  mmdata_ica_p2_exp2%>% mutate(cond = factor(cond, levels = c("DF", "AF")))%>% 
  ggplot(aes(x=targetDur, y=mm_P2, color=cond))+geom_point()+
  geom_errorbar(aes(ymin = mm_P2-se_P2,ymax = mm_P2+se_P2),width = 40, position = position_dodge(width=0.5))+mytheme+
    labs(x = "Target interval (ms)", y = "Mean ampliutde (\u03BCV)") +
    scale_color_manual(values = colors_plot)+ theme(legend.position = c(0.9,0.9), legend.title = element_blank())+ scale_x_continuous(breaks=c(400,800,1200,1600)) + geom_line() + coord_cartesian(ylim = c(-2, 9)) + ggtitle("LPCt")

 m_p2_exp2
```


#2.3- P2 mean amplitude 
```{r}
# test target Dur x condition x mean amplitude
mdata_ica_p2_exp2 <- data_ica_pos_exp2%>% filter( (t>=140 & t<=300)) %>% group_by(targetDur,cond, SubName) %>% dplyr::summarise(n= n(), m_pos = mean(P2),se_pos = sd(P2)/sqrt(n-1))

mdata_ica_p2_exp2_sub= mdata_ica_p2_exp2%>% group_by(cond,SubName) %>% dplyr::summarise(n= n(), mm_pos = mean(m_pos))
mdata_ica_p2_exp2_sub%>% group_by(cond) %>% dplyr::summarise(n= n(), m_pos = mean(mm_pos), se_pos=sd(mm_pos)/sqrt(n-1))
```


###mixed linear model
```{r}
mdata_ica_p2_exp2$targetDur = mdata_ica_p2_exp2$targetDur /1000
```

```{r}
mdata_ica_p2_exp2$cond = as.factor(mdata_ica_p2_exp2$cond)
contrasts(mdata_ica_p2_exp2$cond) = contr.Sum(levels(mdata_ica_p2_exp2$cond))

# fit the linear mixed model
m_lm = lmer(m_pos ~
              cond*targetDur + # fixed effect, covariate and interaction
              (1+targetDur|SubName), # random effect
              data = mdata_ica_p2_exp2)
#output as table
tab_model(m_lm, p.val = 'kr') # Kenward-Roger approximation, alternative 'satterthwaite' approxmiation
```



### Figure 8

```{r}
fig_exp2_pos <- cowplot::plot_grid(plot_1_p2_exp2_df, p2p3_exp2_lat,plot_1_p2_exp2_af,m_p2_exp2, nrow = 2, labels = c("a", "b","c","d"), rel_widths = c(3, 1.3))

ggsave("figures/figure8.png", fig_exp2_pos, width = 8, height = 4.5) 

fig_exp2_pos
```

# CROSS_EXPERIMENT
#plot
```{r}
dat= m_allAverageDat_exp1 %>% filter(targetDur==992|targetDur==1008)
dat2= m_allAverageDat_exp2 %>% filter(targetDur==1000)
datm= rbind(dat,dat2)

 p1000= ggplot(datm, aes(x=t, y=m_CNV, color=as.factor(cond)))+
  scale_color_manual(values=c11)+
  geom_hline(yintercept=0, linetype="dashed")+theme(axis.title.y = element_text(angle=90))+ geom_vline(xintercept = 1000, linetype="dashed")+geom_vline(xintercept = 1300, linetype="dashed")+geom_line(size=0.5)+ 
  scale_y_reverse() + mytheme+ labs(x = "Time (ms)", y = "Amplitude (\u03BCV)", color='') +
  scale_x_continuous(breaks= c(0, 400,800,1200, 1600,2000), labels = c(0, 400,800,1200, 1600,2000)) +coord_cartesian(ylim = c(7.5, -6),xlim = c(-200, 2400), expand = c(0,0))+theme(legend.key.size = unit(0.3, "cm"))+ggtitle("992/1000/1008 ms target intervals")+ theme(legend.position = c(0.9,0.95))+theme(legend.text=element_text(size=8))
 p1000
```


```{r}
dat= m_allAverageDat_exp1 %>% filter(targetDur==400|targetDur==1600)
dat2= m_allAverageDat_exp2 %>% filter(targetDur==400|targetDur==1600)
datm= rbind(dat,dat2)

p1600= ggplot(datm %>% filter(targetDur==1600), aes(x=t, y=m_CNV, color=as.factor(cond)))+
  scale_color_manual(values=c11)+
  scale_y_reverse() + mytheme+ labs(x = "Time (ms)", y = "Amplitude (\u03BCV)", color='') +
  scale_x_continuous(breaks= c(0, 400,800,1200, 1600,2000), labels = c(0, 400,800,1200, 1600,2000)) +coord_cartesian(ylim = c(7.5, -6),xlim = c(-200, 2400), expand = c(0,0))+theme(legend.key.size = unit(0.3, "cm"))+ggtitle("1600 ms target interval")+ theme(legend.position = c(0.95,0.99))+
  geom_hline(yintercept=0, linetype="dashed")+theme(axis.title.y = element_text(angle=90))+ geom_vline(xintercept = 1600, linetype="dashed")+geom_vline(xintercept = 1900, linetype="dashed")+theme(legend.text=element_text(size=8))+geom_line(size=0.5)
p1600
```

```{r}
p400= ggplot(datm %>% filter(targetDur==400), aes(x=t, y=m_CNV, color=as.factor(cond)))+
  scale_color_manual(values=c11)+
  scale_y_reverse() + mytheme+ labs(x = "Time (ms)", y = "Amplitude (\u03BCV)", color='') +
  scale_x_continuous(breaks= c(0, 400,800,1200, 1600,2000), labels = c(0, 400,800,1200, 1600,2000)) +coord_cartesian(ylim = c(7.5, -6),xlim = c(-200, 2400), expand = c(0,0))+theme(legend.key.size = unit(0.3, "cm"))+ggtitle("400 ms target interval")+ theme(legend.position = c(0.9,0.95))+
  geom_hline(yintercept=0, linetype="dashed")+theme(axis.title.y = element_text(angle=90))+ geom_vline(xintercept = 400, linetype="dashed")+geom_vline(xintercept = 700, linetype="dashed")+theme(legend.text=element_text(size=8))+geom_line(size=0.5)
p400
```

#1. CNV slope
# target 400
```{r}
dat_p400= allAverageDat_exp1 %>% filter(targetDur==400)
dat_p400$targetDur= dat_p400$targetDur/1000
dat_p400$t = dat_p400$t /1000
dat2_p400= allAverageDat_exp2 %>% filter(targetDur==400)
dat2_p400$targetDur= dat2_p400$targetDur/1000
dat2_p400$t = dat2_p400$t /1000

dat_p400_ns <- dat_p400 %>% filter(cond == 'NS')
dat_p400_ns %>% filter(t <= 0.65& t >= 0.25) %>% 
  group_by(SubName) %>% nest() %>%
  mutate(model = map(data, cnv_fit)) %>% 
  unnest(model, .drop = TRUE) -> lmdat_p400_ns
lmdat_p400_ns$cond = "NS"
lmdat_p400_ns$exp = "Exp.1"
lmdat_p400_ns$targetDur= 400

dat_p400_ps <- dat_p400 %>% filter(cond == 'PS')
dat_p400_ps %>% filter(t <= 0.65& t >= 0.25) %>% 
  group_by(SubName) %>% nest() %>%
  mutate(model = map(data, cnv_fit)) %>% 
  unnest(model, .drop = TRUE) -> lmdat_p400_ps
lmdat_p400_ps$cond = "PS"
lmdat_p400_ps$exp = "Exp.1"
lmdat_p400_ps$targetDur= 400
```

```{r}
dat2_p400_df <- dat2_p400 %>% filter(cond == 'DF')
dat2_p400_df %>% filter(t <= 0.65& t >= 0.25) %>% 
  group_by(SubName) %>% nest() %>%
  mutate(model = map(data, cnv_fit)) %>% 
  unnest(model, .drop = TRUE) -> lmdat2_p400_df
lmdat2_p400_df$cond = "DF"
lmdat2_p400_df$exp = "Exp.2"
lmdat2_p400_df$targetDur= 400

dat2_p400_af <- dat2_p400 %>% filter(cond == 'AF')
dat2_p400_af %>% filter(t <= 0.65& t >= 0.25) %>% 
  group_by(SubName) %>% nest() %>%
  mutate(model = map(data, cnv_fit)) %>% 
  unnest(model, .drop = TRUE) -> lmdat2_p400_af
lmdat2_p400_af$cond = "AF"
lmdat2_p400_af$exp = "Exp.2"
lmdat2_p400_af$targetDur= 400
```

#target 1600
```{r}
dat_p1600= allAverageDat_exp1 %>% filter(targetDur==1600)
dat_p1600$targetDur= dat_p1600$targetDur/1000
dat_p1600$t = dat_p1600$t /1000
dat2_p1600= allAverageDat_exp2 %>% filter(targetDur==1600)
dat2_p1600$targetDur= dat2_p1600$targetDur/1000
dat2_p1600$t = dat2_p1600$t /1000

dat_p1600_ns <- dat_p1600 %>% filter(cond == 'NS')
dat_p1600_ns %>% filter(t <= 0.65& t >= 0.25) %>% 
  group_by(SubName) %>% nest() %>%
  mutate(model = map(data, cnv_fit)) %>% 
  unnest(model, .drop = TRUE) -> lmdat_p1600_ns
lmdat_p1600_ns$cond = "NS"
lmdat_p1600_ns$exp = "Exp.1"
lmdat_p1600_ns$targetDur=1600

dat_p1600_ps <- dat_p1600 %>% filter(cond == 'PS')
dat_p1600_ps %>% filter(t <= 0.65& t >= 0.25) %>% 
  group_by(SubName) %>% nest() %>%
  mutate(model = map(data, cnv_fit)) %>% 
  unnest(model, .drop = TRUE) -> lmdat_p1600_ps
lmdat_p1600_ps$cond = "PS"
lmdat_p1600_ps$exp = "Exp.1"
lmdat_p1600_ps$targetDur=1600
```

```{r}
dat2_p1600_df <- dat2_p1600 %>% filter(cond == 'DF')
dat2_p1600_df %>% filter(t <= 0.65& t >= 0.25) %>% 
  group_by(SubName) %>% nest() %>%
  mutate(model = map(data, cnv_fit)) %>% 
  unnest(model, .drop = TRUE) -> lmdat2_p1600_df
lmdat2_p1600_df$cond = "DF"
lmdat2_p1600_df$exp = "Exp.2"
lmdat2_p1600_df$targetDur = 1600

dat2_p1600_af <- dat2_p1600 %>% filter(cond == 'AF')
dat2_p1600_af %>% filter(t <= 0.65& t >= 0.25) %>% 
  group_by(SubName) %>% nest() %>%
  mutate(model = map(data, cnv_fit)) %>% 
  unnest(model, .drop = TRUE) -> lmdat2_p1600_af
lmdat2_p1600_af$cond = "AF"
lmdat2_p1600_af$exp = "Exp.2"
lmdat2_p1600_af$targetDur = 1600
```

#target 1000
```{r}
dat_p1000= allAverageDat_exp1 %>% filter(targetDur==1008 | targetDur==992)
dat_p1000$targetDur= dat_p1000$targetDur/1000
dat_p1000$t = dat_p1000$t /1000
dat2_p1000= allAverageDat_exp2 %>% filter(targetDur==1000)
dat2_p1000$targetDur= dat2_p1000$targetDur/1000
dat2_p1000$t = dat2_p1000$t /1000

dat_p1000_ns <- dat_p1000 %>% filter(cond == 'NS')
dat_p1000_ns %>% filter(t <= 0.65& t >= 0.25) %>% 
  group_by(SubName) %>% nest() %>%
  mutate(model = map(data, cnv_fit)) %>% 
  unnest(model, .drop = TRUE) -> lmdat_p1000_ns
lmdat_p1000_ns$cond = "NS"
lmdat_p1000_ns$exp = "Exp.1"
lmdat_p1000_ns$targetDur = 1000

dat_p1000_ps <- dat_p1000 %>% filter(cond == 'PS')
dat_p1000_ps %>% filter(t <= 0.65& t >= 0.25) %>% 
  group_by(SubName) %>% nest() %>%
  mutate(model = map(data, cnv_fit)) %>% 
  unnest(model, .drop = TRUE) -> lmdat_p1000_ps
lmdat_p1000_ps$cond = "PS"
lmdat_p1000_ps$exp = "Exp.1"
lmdat_p1000_ps$targetDur = 1000
```

```{r}
dat2_p1000_df <- dat2_p1000 %>% filter(cond == 'DF')
dat2_p1000_df %>% filter(t <= 0.65& t >= 0.25) %>% 
  group_by(SubName) %>% nest() %>%
  mutate(model = map(data, cnv_fit)) %>% 
  unnest(model, .drop = TRUE) -> lmdat2_p1000_df
lmdat2_p1000_df$cond = "DF"
lmdat2_p1000_df$exp = "Exp.2"
lmdat2_p1000_df$targetDur = 1000

dat2_p1000_af <- dat2_p1000 %>% filter(cond == 'AF')
dat2_p1000_af %>% filter(t <= 0.65& t >= 0.25) %>% 
  group_by(SubName) %>% nest() %>%
  mutate(model = map(data, cnv_fit)) %>% 
  unnest(model, .drop = TRUE) -> lmdat2_p1000_af
lmdat2_p1000_af$cond = "AF"
lmdat2_p1000_af$exp = "Exp.2"
lmdat2_p1000_af$targetDur = 1000
```

```{r}
lm= rbind(lmdat2_p1000_df, lmdat2_p1000_af,lmdat_p1000_ps,lmdat_p1000_ns,
          lmdat2_p400_df, lmdat2_p400_af,lmdat_p400_ps,lmdat_p400_ns,
          lmdat2_p1600_df, lmdat2_p1600_af,lmdat_p1600_ps,lmdat_p1600_ns)
lm$Sub = lm$SubName
lm%>% group_by(exp,cond) %>% dplyr::summarise(m_slope= mean(slope),n=n(), se_slope=sd(slope)/sqrt(n-1))
```

```{r}
lm$cond[lm$cond =='AF']='long'
lm$cond[lm$cond =='DF']='short'
lm$cond[lm$cond =='PS']='short'
lm$cond[lm$cond =='NS']='long'
```

#lmm

```{r}
lm$cond = as.factor(lm$cond)
contrasts(lm$cond) = contr.Sum(levels(lm$cond))
```

```{r}
# fit the linear mixed model
lmm_comm = lmer(slope ~exp*cond + # fixed effect, covariate and interaction
              (targetDur|SubName), # random effect
              data = lm)
#output as table
tab_model(lmm_comm, p.val = 'kr') # Kenward-Roger approximation, alternative 'satterthwaite' approxmiation
```


# 2. CNV end
```{r}
m_turn_exp1= turn_exp1%>% group_by(targetDur, cond) %>% dplyr::summarise(n=n(),point= mean(time), se_point= sd(time)/sqrt(n-1))
m_turn_exp1$cond[m_turn_exp1$cond=="NS"]= "long"
m_turn_exp1$cond[m_turn_exp1$cond=="PS"]= "short"
m_turn_exp1$targetDur[m_turn_exp1$targetDur==1008]= 1000
m_turn_exp1$targetDur[m_turn_exp1$targetDur==992]= 1000
m_turn_exp1$Exp= "Exp.1"

m_turn_exp2= turn_exp2%>% group_by(targetDur, cond) %>% dplyr::summarise(n=n(),point= mean(time), se_point= sd(time)/sqrt(n-1))
m_turn_exp2$cond[m_turn_exp2$cond=="AF"]= "long"
m_turn_exp2$cond[m_turn_exp2$cond=="DF"]= "short"
m_turn_exp2$Exp= "Exp.2"

dat_turn = rbind(m_turn_exp1, m_turn_exp2)
dat_turn %>% group_by(Exp, targetDur) %>% dplyr::summarise(n=n(),m_point=mean(point), se_point= sd(point)/sqrt(n-1))
```

#t-test

```{r}
turn_exp1$cond[turn_exp1$cond=="NS"]= "long"
turn_exp1$cond[turn_exp1$cond=="PS"]= "short"
turn_exp1$targetDur[turn_exp1$targetDur==1008]= 1000
turn_exp1$targetDur[turn_exp1$targetDur==992]= 1000
turn_exp1$Exp= "Exp.1"
turn_exp2$cond[turn_exp2$cond=="AF"]= "long"
turn_exp2$cond[turn_exp2$cond=="DF"]= "short"
turn_exp2$Exp= "Exp.2"

dat_turn_all = rbind(turn_exp1, turn_exp2)
t4= turn_exp1%>%filter(targetDur==1600)
t42=turn_exp2%>%filter(targetDur==1600)

t.test(t4$time, t42$time)
```

```{r}
t3= turn_exp1%>%filter(targetDur==1000)
t32=turn_exp2%>%filter(targetDur==1000)

t.test(t3$time, t32$time)
```
```{r}
t2= turn_exp1%>%filter(targetDur==400)
t22=turn_exp2%>%filter(targetDur==400)

t.test(t2$time, t22$time)
```

#plot
```{r}
col1=  c("#000333", "#3399CC", "#FF3333")

mm_turn_exp1= dat_turn%>% group_by(targetDur,Exp) %>% dplyr::summarise(n=n(),mpoint= mean(point), se_point= sd(point)/sqrt(n-1))
mm_turn_exp1$Exp= as.factor(mm_turn_exp1$Exp)
mm_turn_exp1$targetDur= as.factor(mm_turn_exp1$targetDur)
p_cross= ggplot(data=mm_turn_exp1, aes(x=Exp, y=mpoint, color=targetDur,fill=targetDur))+ geom_hline(yintercept = 400,linetype="dashed")+ geom_hline(yintercept = 1000,linetype="dashed")+ geom_hline(yintercept = 1600,linetype="dashed") +geom_bar(stat="identity",position=position_dodge(),width=0.8,fill="white")+
  mytheme+ scale_color_manual(values=col1)+ theme(legend.position = "top") + labs(x="", y="CNV end (ms)",color="",fill="") +  geom_errorbar(aes(ymin=mpoint-se_point, ymax=mpoint+se_point, color=targetDur), width=.2,position=position_dodge(0.8))+theme(legend.text=element_text(size=7))+theme(legend.key.size = unit(0.2, "cm")) +theme(axis.text.x = element_text(size=8)) +theme(axis.text.y = element_text(size=8),axis.title.y = element_text(size=8))
p_cross
```


#barplot

```{r}
lm= rbind(lmdat2_p1000_df, lmdat2_p1000_af,lmdat_p1000_ps,lmdat_p1000_ns,
          lmdat2_p400_df, lmdat2_p400_af,lmdat_p400_ps,lmdat_p400_ns,
          lmdat2_p1600_df, lmdat2_p1600_af,lmdat_p1600_ps,lmdat_p1600_ns)
lm$Sub = lm$SubName
lm_m = lm%>% group_by(cond, exp, SubName) %>% dplyr::summarise(mslope=mean(slope))
lm_mm=lm_m%>% group_by(cond, exp) %>% dplyr::summarise(m_slope= mean(mslope), n=n(), se_slope=sd(mslope)/sqrt(n-1))

p_slope= ggplot(data=lm_mm, aes(x=exp, y=m_slope, color=cond,fill=cond)) +geom_bar(stat="identity",position=position_dodge(), width=0.5, fill="white")+ scale_y_reverse() +mytheme+ scale_color_manual(values=c11)+ theme(legend.position = "top") + labs(x="", y="CNV slope (s/\u03BCV)",color="",fill="") + geom_errorbar(aes(ymin = m_slope-se_slope,ymax = m_slope+ se_slope),width = 0.1, position = position_dodge(width=0.5)) +theme(legend.text=element_text(size=7))+ coord_cartesian(ylim=c(0,-23))+theme(legend.key.size = unit(0.2, "cm"))+theme(axis.text.x = element_text(size=8)) +theme(axis.text.y = element_text(size=8),axis.title.y = element_text(size=8))
p_slope
```


### Figure 9

```{r}
lastrow= cowplot::plot_grid(p1600,p_slope,p_cross,labels=c("c","d","e"),rel_widths = c(1,0.5,0.5),nrow=1)
firstrow <- cowplot::plot_grid(p400,p1000,nrow = 1, labels = c("a", "b"), rel_widths = c(1,1))
cross_plot =cowplot::plot_grid(firstrow, lastrow, rel_heights = c(1,1),nrow=2)

ggsave("figures/figure9.png", cross_plot, width = 8, height = 4.9) 

cross_plot
```