EEG_temporal_reproduction/derivatives/EEG_analysis.Rmd

---
title: "Blocked vs. Interleaved: How range contexts modulate time perception and its EEG signatures"
author: "Cemre Baykan, Xiuna Zhu, Artyom Zinchenko, Zhuanghua Shi"
date: "August 2023"
---

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

### 1. Perception
```{r}
# plot
p1t_Perception_all <- ggplot(m_validAverageDat_short %>% filter(phase=="Perception") ,aes(x = t,y = CNV_m, colour = as.factor(targetDur), linetype=cond)) +geom_vline(xintercept = 300, linetype="dashed")+geom_vline(xintercept = 600, linetype="dashed")+
  geom_line(na.rm = FALSE, size=0.5) + labs(x = "Time (s)", y = "Amplitude(\u03BCV)", color='', linetype='') + 
  geom_vline(xintercept = 0, size=0.2, color="grey") +
  geom_hline(yintercept = 0, size=0.2,color="grey") +
  coord_cartesian(ylim = c(8, -8), xlim = c(-200, 3600), expand = c(0,0))+theme_new2+ theme(legend.position = "right" )+ scale_x_continuous(breaks=c(0,400,1200,2400,3600), labels=c(0,0.4,1.2,2.4,3.6))+scale_color_manual(values = colors3,labels=c("0.4","0.566","0.8","1.2", "1.697","2.4"))
p1t_Perception_all
```

###1.1 CNV buildup rate
```{r}
allAverageDat1 = allAverageDat%>% filter( (t>= 300 & t <= 600)) %>% dplyr::group_by(cond,group, phase, SubName) %>% dplyr::summarise(n= n(), m_cnv = mean(CNV))

allAverageDat1_sub = allAverageDat1 %>% dplyr::group_by(cond,SubName, phase) %>% dplyr::summarise(n= n(), mm_cnv = mean(m_cnv))

allAverageDat1_sub %>% dplyr::group_by(cond, phase) %>% dplyr::summarise(n= n(), mmm_cnv = mean(mm_cnv), se_cnv= sd(mm_cnv)/sqrt(n-1))
```

```{r}
mallAverageDat1 <- allAverageDat1 %>% dplyr::group_by(phase,cond,group) %>% dplyr::summarise(mm_cnv = mean(m_cnv), se_cnv = sd(m_cnv)/sqrt(n-1), n= n())
mallAverageDat1=unique(mallAverageDat1)
```


```{r}
mallAverageDat1_per= mallAverageDat1%>% filter(phase=="Perception")
targ_labeller2=c("Lower","Upper")
p_cnv1_target= ggplot(data=mallAverageDat1_per%>% mutate(cond = factor(cond, levels = c("BR", "IR")))%>% mutate(group = factor(group, levels = c("Lower", "Upper"))), aes(x=group, y=mm_cnv, color=group,fill=group,linetype=cond))+ geom_hline(yintercept = -0.5, linetype="dashed")  +geom_hline(yintercept = -1, linetype="dashed")+ geom_hline(yintercept = 0, linetype="dashed") +geom_bar(stat="identity",position=position_dodge(), width=0.6, fill="white", alpha=0.9)+ theme_new + scale_color_manual(values=col11)+ theme(legend.position ="bottom") + labs(x="", y="Amplitude (\u03BCV)",color="",linetype="") + geom_errorbar(aes(ymin = mm_cnv-se_cnv,ymax = mm_cnv+ se_cnv),width = 0.1, position = position_dodge(width=0.6))+ scale_y_reverse()+  ggtitle("CNV buildup")+guides(colour = "none")+ coord_cartesian(ylim=c(0,-3)) +theme(legend.text=element_text(size=9))+theme(legend.key.size = unit(0.5, "cm"))+theme(legend.box.spacing = unit(-10, "pt"))

p_cnv1_target
```


# anova
```{r}
mallAverageDat1_sub <- allAverageDat1 %>% group_by(phase,cond,SubName,group) %>% dplyr::summarise(mm_cnv = mean(m_cnv))

ezANOVA(data = allAverageDat1%>% filter(phase=="Perception"), dv= m_cnv, wid=SubName, within=.(cond, group))
```
### 1.2 CNV mean ampl.

```{r}
# calculate average CNV until offsets
cnv_build2 <- data.frame()
targetList <- c(400,566,800,1200,1697,2400)
for (target in targetList){
  m_cnv_1 <- allAverageDat %>%filter(t>= 300 & t <= (target+300))%>% dplyr::group_by(cond,group, phase, SubName, targetDur) %>% dplyr::summarise(n= n(), m_CNV = mean(CNV))
  m_cnv_1= m_cnv_1%>% filter(targetDur==target)
  cnv_build2 <- rbind2(cnv_build2,m_cnv_1, by=c("SubName", "group","phase","targetDur","cond")) 
}
```

```{r}
mallAverageDat2 <- cnv_build2 %>% dplyr::group_by(phase,cond, group,targetDur) %>% dplyr::summarise(mm_cnv = mean(m_CNV), se_cnv = sd(m_CNV)/sqrt(n-1), n= n())
mallAverageDat2= unique(mallAverageDat2)
tallAverageDat2 <- mallAverageDat2 %>% dplyr::group_by(phase,cond, group) %>% dplyr::summarise(mmm_cnv = mean(mm_cnv), se_cnv = sd(mm_cnv)/sqrt(n-1), n= n())
tallAverageDat2= unique(tallAverageDat2)
```  


```{r}
targ_labeller= c('400'='0.4', '566'='0.566','800'='0.8','1200'='1.2','1697'='1.697','2400'='2.4')

mallAverageDat2$targetDur= as.factor(mallAverageDat2$targetDur)
p_cnv2_target= ggplot(data=mallAverageDat2%>% filter(phase=="Perception")%>% mutate(cond = factor(cond, levels = c("BR", "IR")))%>% mutate(group = factor(group, levels = c("Lower", "Upper"))), aes(x=cond, y=mm_cnv, color=group,fill=group, linetype=cond))+ geom_hline(yintercept = -1, linetype="dashed")+geom_hline(yintercept = -2, linetype="dashed")+geom_hline(yintercept = -3, linetype="dashed")+ geom_hline(yintercept = 0, linetype="dashed") +geom_bar(stat="identity",position=position_dodge(), width=0.6, fill="white", alpha=0.9)+ theme_new + theme(legend.position = "bottom") + labs(x="", y="Amplitude (\u03BCV)",color="", linetype="") + geom_errorbar(aes(ymin = mm_cnv-se_cnv,ymax = mm_cnv+ se_cnv),width = 0.1, position = position_dodge(width=0.6)) + scale_y_reverse() + scale_color_manual(values=col11) + ggtitle("Mean CNV")+ scale_x_discrete(labels=c(0.4,0.566,0.8, 1.2,1.697,2.4))+facet_grid(~targetDur,labeller=as_labeller(targ_labeller))+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x = element_blank())+theme(legend.box.spacing = unit(3, "pt"))

p_cnv2_target
```

# lmm
```{r}
cnv_build2_p= cnv_build2 %>% filter(phase=="Perception")
cnv_build2_p= cnv_build2_p%>% dplyr::mutate(cond= factor(cond,levels=c("IR","BR")))
cnv_build2_p= cnv_build2_p%>% dplyr::mutate(group= factor(group,levels=c("Lower","Upper")))
cnv_build2_p$targetDur_n= cnv_build2_p$targetDur-1200
cnv_build2_p$targetDur_n= cnv_build2_p$targetDur_n/1000
cnv_build2$cond = as.factor(cnv_build2$cond)
cnv_build2$group = as.factor(cnv_build2$group)
contrasts(cnv_build2$cond) = contr.Sum(levels(cnv_build2$cond))
contrasts(cnv_build2$group) = contr.Sum(levels(cnv_build2$group))
```

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

## 1.3 CNV width
```{r}
mcnv_width_per <- cnv_width_per %>% group_by(cond, targetDur) %>% dplyr::summarise(n= n(),mlat = mean(time), se_lat = sd(time)/sqrt(n-1))
```
#lmm

```{r}
cnv_width_per$cond = as.factor(cnv_width_per$cond)
contrasts(cnv_width_per$cond) = contr.Sum(levels(cnv_width_per$cond))
cnv_width_per$targetDur_n = cnv_width_per$targetDur - 1200
```

```{r}
scnv_width_per <- cnv_width_per %>% group_by(cond, targetDur,targetDur_n, SubName,group) %>% dplyr::summarise(n= n(),time = mean(time))
```

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

```{r}
m_scnv_width_per = scnv_width_per%>% group_by(cond, SubName) %>% dplyr::summarise(n= n(),time = mean(time))
m_scnv_width_per %>% group_by(cond) %>% dplyr::summarise(time = mean(time))
```

######figure 3
```{r}
frow <- plot_grid(p_cnv1_target,p_cnv2_target, nrow=1, rel_widths  = c(1, 3), labels=c("c","d"))

srow <- plot_grid(NULL,p1t_Perception_all,frow, nrow=3, rel_heights  = c(0.5,1,0.8), labels=c("a","b",""))

ggsave("figures/figure3.png", srow, width = 8.5, height = 6.5) 
srow
```

## 2. Offset-P2

```{r}
mallAverageDat_p2=allAverageDat_p2%>% group_by(phase, group, cond, targetDur,t) %>% dplyr::summarise(mP2 = mean(P2), sd_P2 = sd(P2),n= n())
```


```{r}
#mix condition
p1t_Perception_p2 <- ggplot(mallAverageDat_p2 %>% filter(phase=="Perception") ,aes(x = t,y = mP2, colour = as.factor(targetDur), linetype=cond))+
  geom_line(na.rm = FALSE) + labs(x = "Time (s)", y = "Amplitude(\u03BCV)", color='', linetype='') + 
  geom_vline(xintercept = 0, size=0.2, color="grey") +
  geom_hline(yintercept = 0, size=0.2,color="grey") +
  scale_x_continuous(breaks= c(0, 200, 600,1000), labels=c("0","0.2","0.6","1.0")) +theme_new2+ theme(legend.key.size = unit(0.5, "cm"))+scale_color_manual(values = colors3,labels=c("0.4","0.566","0.8","1.2", "1.697","2.4"))+
  scale_y_reverse()+coord_cartesian(ylim = c(8, -5), xlim = c(-100, 1000), expand = c(0,0))+theme(legend.position = "right") 
  
p1t_Perception_p2
```


### 2.1 P2 mean ampl.

```{r}
m_validAverageDat_p2 = allAverageDat_p2%>% filter( (t>= 200 & t <= 300)) %>% group_by(cond,group,targetDur, phase, SubName) %>% dplyr::summarise(n= n(), m_P2 = mean(P2),se_P2 = sd(P2)/sqrt(n-1))

m_validAverageDat_p2_percep= m_validAverageDat_p2 %>% filter(phase=="Perception")
m_validAverageDat_p2_reprod= m_validAverageDat_p2 %>% filter(phase=="Reproduction")
```

```{r}
mm_validAverageDat_p2 <- m_validAverageDat_p2 %>% group_by(phase,cond, group,targetDur) %>% dplyr::summarise(mm_P2 = mean(m_P2), se_P2 = sd(m_P2)/sqrt(n-1), n= n())

mm_validAverageDat_p2_percep= mm_validAverageDat_p2 %>% filter(phase=="Perception")
mm_validAverageDat_p2_reprod= mm_validAverageDat_p2 %>% filter(phase=="Reproduction")
```  
  

```{r}
targ_labeller2= c('0.4'='0.4', '0.566'='0.566','0.8'='0.8','1.2'='1.2','1.697'='1.697','2.4'='2.4')

p_cnv2_target_p2= ggplot(data=mm_validAverageDat_p2_percep%>% filter(phase=="Perception")%>% mutate(cond = factor(cond, levels = c("BR", "IR")))%>% mutate(group = factor(group, levels = c("Lower", "Upper"))), aes(x=cond, y=mm_P2, color=group,fill=group,linetype=cond)) +geom_bar(stat="identity",position=position_dodge(), width=0.6, fill="white", alpha=0.9)+ theme_new2+ scale_color_manual(values=col11)+ theme(legend.position = "bottom") + labs(x="", y="Amplitude (\u03BCV)",color="",linetype="") + geom_errorbar(aes(ymin = mm_P2-se_P2,ymax = mm_P2+ se_P2),width = 0.1, position = position_dodge(width=0.6)) + ggtitle("Mean P2") +facet_grid(~targetDur,labeller=as_labeller(targ_labeller2))+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())+theme(legend.key.size = unit(0.3, "cm")) 

p_cnv2_target_p2
```
#lmm
```{r}
m_validAverageDat_p2_percep$cond = as.factor(m_validAverageDat_p2_percep$cond)
m_validAverageDat_p2_percep$group = as.factor(m_validAverageDat_p2_percep$group)
contrasts(m_validAverageDat_p2_percep$cond) = contr.Sum(levels(m_validAverageDat_p2_percep$cond))
contrasts(m_validAverageDat_p2_percep$group) = contr.Sum(levels(m_validAverageDat_p2_percep$group))
m_validAverageDat_p2_percep$targetDur_n = m_validAverageDat_p2_percep$targetDur - 1.2
```

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

```{r}
mm_validAverageDat_p2_percep  %>% group_by(cond) %>% dplyr::summarise(mm_P2 = mean(mm_P2),n=n())
```
###### figure 4
```{r}
frow_31 <- plot_grid(NULL,p1t_Perception_p2,nrow=2,ncol = 1, rel_heights  = c(0.35,1), labels=c("a",""))

frow_3 <- plot_grid(frow_31,p_cnv2_target_p2,nrow=1,ncol = 2, rel_widths  = c(3.2, 1.7), labels=c("","b"))

ggsave("figures/figure4.png", frow_3, width = 8.5, height = 2.7) 
frow_3
```

# 3. LPCt
```{r}
m_allAverageDat_resp <- allAverageDat_resp %>% group_by(phase, group, cond, targetDur,t) %>% dplyr::summarise(m_p3 = mean(P3), sd_p3 = sd(P3),n= n())
```


```{r}
plt_resp <- ggplot(m_allAverageDat_resp %>% filter(phase=="Perception") ,aes(x = t,y = m_p3, colour = as.factor(targetDur), linetype=cond)) +
  geom_line(na.rm = FALSE, size=0.5) + labs(x = "Time (s)", y = "Amplitude(\u03BCV)", color='', linetype='') +
  geom_vline(xintercept = 0, size=0.2, color="grey") +
  geom_hline(yintercept = 0, size=0.2,color="grey") +theme_new2+ theme(legend.position = "right" )+theme(legend.key.size = unit(0.5, "cm"))+ scale_x_continuous(breaks=c(-1300,-800,-300,0,300), labels=c(-1.3,-0.8,-0.3,0,0.3)) + scale_y_reverse() +coord_cartesian(ylim = c(8, -8), expand = c(0,0)) + scale_color_manual(values = colors3,labels=c("0.4","0.566","0.8","1.2", "1.697","2.4"))+ coord_cartesian(xlim=c(-1000,500), ylim=c(7,-5))
  
plt_resp
```

```{r}
allAverageDat_resp$targetDur = allAverageDat_resp$targetDur* 1000
```

```{r}
# calculate average amplitude in pre-reproduction
cnv_build3 <- data.frame()
targetList <- c(400,566,800,1200,1697,2400)
for (target in targetList){
  m_cnv_1 <- allAverageDat_resp %>%filter(t>= (-800) & t <= (-600))%>% group_by(cond,group, phase, SubName, targetDur) %>% dplyr::summarise(n= n(), m_P3 = mean(P3),se_P3 = sd(P3)/sqrt(n-1))
  m_cnv_1= m_cnv_1%>% filter(targetDur==target)
  cnv_build3 <- rbind2(cnv_build3,m_cnv_1, by=c("SubName", "group","phase","targetDur","cond")) 
}
```

```{r}
mallAverageDat3 <- cnv_build3 %>% group_by(phase,cond, group,targetDur) %>% dplyr::summarise(mm_p3 = mean(m_P3), se_p3 = sd(m_P3)/sqrt(n-1), n= n())
mallAverageDat3= unique(mallAverageDat3)

tallAverageDat3 <- mallAverageDat3 %>% group_by(phase,cond, group) %>% dplyr::summarise(mmm_p3 = mean(mm_p3), se_p3 = sd(mm_p3)/sqrt(n-1), n= n())
tallAverageDat3= unique(tallAverageDat3)
```  


```{r}
targ_labeller= c('400'='0.4', '566'='0.566','800'='0.8','1200'='1.2','1697'='1.697','2400'='2.4')

p_cnv3_target= ggplot(data=mallAverageDat3%>% filter(phase=="Perception")%>% mutate(cond = factor(cond, levels = c("BR", "IR")))%>% mutate(group = factor(group, levels = c("Lower", "Upper"))), aes(x=cond, y=mm_p3, color=group,fill=group,linetype=cond)) +geom_bar(stat="identity",position=position_dodge(), width=0.6, fill="white", alpha=0.9)+ theme_new2+ scale_color_manual(values=col11)+ theme(legend.position = "bottom") + labs(x="", y="Amplitude (\u03BCV)",color="",linetype="") + geom_errorbar(aes(ymin = mm_p3-se_p3,ymax = mm_p3+ se_p3),width = 0.1, position = position_dodge(width=0.6)) + ggtitle("Mean LPCt") +facet_grid(~targetDur,labeller=as_labeller(targ_labeller))+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())+theme(legend.key.size = unit(0.3, "cm")) 

p_cnv3_target
```

#lmm

```{r}
cnv_build3_p = cnv_build3%>%filter(phase=='Perception')
cnv_build3_p$cond = as.factor(cnv_build3_p$cond)
contrasts(cnv_build3_p$cond) = contr.Sum(levels(cnv_build3_p$cond))
cnv_build3_p$targetDur=cnv_build3_p$targetDur/1000
cnv_build3_p$targetDur_n = cnv_build3_p$targetDur - 1.2
cnv_build3_p$group = as.factor(cnv_build3_p$group)
```

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


```{r}
tallAverageDat3%>% filter(phase=="Perception")%>% group_by(cond)%>% dplyr::summarise(m_p3= mean(mmm_p3))
```

```{r}
tallAverageDat3%>% filter(phase=="Perception")%>% group_by(group,cond)%>% dplyr::summarise(m_p3= mean(mmm_p3))
```

######figure 5

```{r}
frow_21 <- plot_grid(NULL,plt_resp,nrow=2,ncol = 1, rel_heights  = c(0.35, 1), labels=c("a",""))
                     
frow_2 <- plot_grid(frow_21,p_cnv3_target,nrow=1,ncol = 2, rel_widths  = c(3.2, 1.8), labels=c("","b"))

ggsave("figures/figure5.png", frow_2, width = 8.5, height = 2.7) 
frow_2
```

### 4. Reproduction

```{r}
m_validAverageDat <- allAverageDat %>% group_by(phase, group, cond, targetDur,t) %>% dplyr::summarise(CNV_m = mean(CNV), CNV_sd = sd(CNV),n= n())

m_validAverageDat_split <- split(m_validAverageDat, f = list(m_validAverageDat$phase,m_validAverageDat$cond))
#full range condition
p1t_Reproduction_all <- ggplot(m_validAverageDat_short %>% filter(phase=="Reproduction") ,aes(x = t,y = CNV_m, colour = as.factor(targetDur), linetype=cond)) +
  geom_line(na.rm = FALSE, size=0.5) + labs(x = "Time (s)", y = "Amplitude(\u03BCV)", color='', linetype='') +geom_vline(xintercept = 300, linetype="dashed")+geom_vline(xintercept = 600, linetype="dashed")+
  geom_vline(xintercept = 0, size=0.2, color="grey") +
  geom_hline(yintercept = 0, size=0.2,color="grey") +
  scale_x_continuous(breaks= c(0, 400, 1400, 2400)) + 
  coord_cartesian(ylim = c(8, -8), xlim = c(-200, 3600), expand = c(0,0))+theme_new2+  scale_x_continuous(breaks=c(0,400,1200,2400,3600), labels=c(0,0.4,1.2,2.4, 3.6))+ scale_color_manual(values = colors3, labels=c("0.4", "0.566","0.8", "1.2","1.697","2.4"))
  
p1t_Reproduction_all
```
## 4.1 CNV buildup rate

```{r}
p_cnv1_target_p= ggplot(data=mallAverageDat1%>% filter(phase=="Reproduction")%>% mutate(cond = factor(cond, levels = c("BR", "IR")))%>% mutate(group = factor(group, levels = c("Lower", "Upper"))), aes(x=group, y=mm_cnv, color=group,fill=group,linetype=cond))+ geom_hline(yintercept = -2, linetype="dashed") +geom_hline(yintercept = -1, linetype="dashed")+ geom_hline(yintercept = 0, linetype="dashed") +geom_bar(stat="identity",position=position_dodge(), width=0.6, fill="white", alpha=0.9)+ theme_new + scale_color_manual(values=col11)+ theme(legend.position = "top") + labs(x="", y="Amplitude (\u03BCV)",color="", linetype='') + geom_errorbar(aes(ymin = mm_cnv-se_cnv,ymax = mm_cnv+ se_cnv),width = 0.1, position = position_dodge(width=0.6)) + scale_y_reverse() +theme(legend.key.size = unit(0.3, "cm"))+ coord_cartesian(ylim=c(0,-3.5)) + ggtitle("CNV buildup") +guides(color = "none")+ guides(linetype=guide_legend(nrow=2,byrow=TRUE))+ theme(legend.position = c(0.75,0.9))

p_cnv1_target_p
```

```{r}
allAverageDat1_sub2 = allAverageDat1 %>% dplyr::group_by(group,SubName, phase) %>% dplyr::summarise(n= n(), mm_cnv = mean(m_cnv))

allAverageDat1_sub2 %>% dplyr::group_by(group, phase) %>% dplyr::summarise(n= n(), mmm_cnv = mean(mm_cnv), se_cnv= sd(mm_cnv)/sqrt(n-1))
```

#anova

```{r}
ezANOVA(data = allAverageDat1%>% filter(phase=="Reproduction"), dv= m_cnv, wid=SubName, within=.(cond, group))
```
###4.2 CNV mean ampl.

#lmm
```{r}
cnv_build2_r= cnv_build2%>% filter(phase=="Reproduction")
cnv_build2_r$targetDur=cnv_build2_r$targetDur/1000
cnv_build2_r$targetDur_n = cnv_build2_r$targetDur - 1.2
```

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

### 4.3 CNV width

```{r}
mcnv_width_rep <- cnv_width_rep %>% group_by(cond, targetDur) %>% dplyr::summarise(n= n(),mlat = mean(time), se_lat = sd(time)/sqrt(n-1))
```
```{r}
mcnv_width_rep %>% group_by(cond) %>% dplyr::summarise(n= n(),time = mean(mlat))
```

```{r}
cnv_width_rep$cond = as.factor(cnv_width_rep$cond)
contrasts(cnv_width_rep$cond) = contr.Sum(levels(cnv_width_rep$cond))
cnv_width_rep$targetDur=cnv_width_rep$targetDur*1000
cnv_width_rep$targetDur_n = cnv_width_rep$targetDur - 1200
```

```{r}
scnv_width_rep <- cnv_width_rep %>% group_by(cond, targetDur,targetDur_n, SubName,group) %>% dplyr::summarise(n= n(),time = mean(time))
scnv_width_rep$targetDur_n = scnv_width_rep$targetDur- 1200
```
#lmm
```{r}
scnv_width_rep= scnv_width_rep%>% dplyr::mutate(cond= factor(cond,levels=c("IR","BR")))
scnv_width_rep= scnv_width_rep%>% dplyr::mutate(group= factor(group,levels=c("Lower","Upper")))
# fit the linear mixed model
fit_m2_p = lmer(time ~
              cond*targetDur_n +group+ # fixed effect, covariate and interaction
              (targetDur_n|SubName), # random effect
              data = scnv_width_rep)
#output as table
tab_model(fit_m2_p, p.val = 'kr') # Kenward-Roger approximation, alternative 'satterthwaite' approxmiation
```

# figure 6
```{r}
frow_r <- plot_grid(NULL,p1t_Reproduction_all,nrow=2,ncol = 1, rel_heights  = c(0.35,1), labels=c("a",""))

frow_r2 <- plot_grid(frow_r,p_cnv1_target_p,nrow=1,ncol = 2, rel_widths  = c(4.5,1), labels=c("","b"))

ggsave("figures/figure6.png", frow_r2, width = 8.5, height = 3.2) 
frow_r2
```