Upload files to 'derivatives'

derivatives/EEG_analysis.Rmd

@@ -1,9 +1,7 @@
 ---
-title: Exploring range effects in the temporal reproduction using EEG signatures of
-  time
+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"
-output: pdf_document
 ---
 
 ```{r setup, include=FALSE}
@@ -24,26 +22,27 @@ p1t_Perception_all
 
 ###1.1 CNV buildup rate
 ```{r}
-allAverageDat1 = allAverageDat%>% filter( (t>= 300 & t <= 600)) %>% group_by(cond,group,targetDur, phase, SubName) %>% dplyr::summarise(n= n(), m_cnv = mean(CNV))
+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 %>% group_by(phase,cond,targetDur,group) %>% dplyr::summarise(mm_cnv = mean(m_cnv), se_cnv = sd(m_cnv)/sqrt(n-1), n= n())
-
-tallAverageDat1 <- mallAverageDat1 %>% group_by(phase,cond, group) %>% dplyr::summarise(mmm_cnv = mean(mm_cnv), se_cnv = sd(mm_cnv)/sqrt(n-1), n= n())
+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=tallAverageDat1%>% filter(phase=="Perception")%>% mutate(cond = factor(cond, levels = c("BR", "IR")))%>% mutate(group = factor(group, levels = c("Lower", "Upper"))), aes(x=group, y=mmm_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 = mmm_cnv-se_cnv,ymax = mmm_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= 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
 ```
 
-```{r}
-tallAverageDat1%>% filter(phase=="Perception")%>% group_by(cond)%>% dplyr::summarise(mcnv= mean(mmm_cnv))
-```
 
 # anova
 ```{r}
@@ -86,18 +85,22 @@ 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))
-cnv_build2_p= cnv_build2%>% filter(phase=="Perception")
-#cnv_build2_p$targetDur=cnv_build2_p$targetDur/1000 
-cnv_build2_p$Dur = cnv_build2_p$targetDur - 1200
+contrasts(cnv_build2$group) = contr.Sum(levels(cnv_build2$group))
 ```
 
 ```{r}
 # fit the linear mixed model
 fit_m2 = lmer(m_CNV ~
-              cond*Dur +group+ # fixed effect, covariate and interaction
-              (Dur|SubName), # random effect
+              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
@@ -112,27 +115,26 @@ mcnv_width_per <- cnv_width_per %>% group_by(cond, targetDur) %>% dplyr::summari
 ```{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=cnv_width_per$targetDur/1000
-cnv_width_per$Dur = cnv_width_per$targetDur - 1200
+cnv_width_per$targetDur_n = cnv_width_per$targetDur - 1200
 ```
 
 ```{r}
-scnv_width_per <- cnv_width_per %>% group_by(cond, targetDur,Dur, SubName,group) %>% dplyr::summarise(n= n(),time = mean(time))
-#scnv_width_per$time= scnv_width_per$time/1000
+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*Dur+group + # fixed effect, covariate and interaction
-              (Dur|SubName), # random effect
+              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}
-scnv_width_per %>% group_by(cond) %>% dplyr::summarise(n= n(),time = mean(time))
+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
@@ -149,7 +151,6 @@ srow
 
 ```{r}
 mallAverageDat_p2=allAverageDat_p2%>% group_by(phase, group, cond, targetDur,t) %>% dplyr::summarise(mP2 = mean(P2), sd_P2 = sd(P2),n= n())
-
 ```
 
 
@@ -196,23 +197,24 @@ 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))
-m_validAverageDat_p2_percep$targetDur=m_validAverageDat_p2_percep$targetDur*1000
-m_validAverageDat_p2_percep$Dur = m_validAverageDat_p2_percep$targetDur - 1200
+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*Dur+group + # fixed effect, covariate and interaction
-              (Dur|SubName), # random effect
+              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))
+mm_validAverageDat_p2_percep  %>% group_by(cond) %>% dplyr::summarise(mm_P2 = mean(mm_P2),n=n())
 ```
 ###### figure 4
 ```{r}
@@ -274,25 +276,22 @@ p_cnv3_target= ggplot(data=mallAverageDat3%>% filter(phase=="Perception")%>% mut
 p_cnv3_target
 ```
 
-```{r}
-cnv_build3_p = cnv_build3%>%filter(phase=='Perception')
-```
-
 #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$Dur = cnv_build3_p$targetDur - 1200
+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*Dur+group+   # fixed effect, covariate and interaction
-              (Dur|SubName), # random effect
+              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
@@ -337,12 +336,15 @@ p1t_Reproduction_all
 ## 4.1 CNV buildup rate
 
 ```{r}
-p_cnv1_target_p= ggplot(data=tallAverageDat1%>% filter(phase=="Reproduction")%>% mutate(cond = factor(cond, levels = c("BR", "IR")))%>% mutate(group = factor(group, levels = c("Lower", "Upper"))), aes(x=group, y=mmm_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 = mmm_cnv-se_cnv,ymax = mmm_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= 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}
-tallAverageDat1%>% filter(phase=="Reproduction")%>% group_by(group)%>% dplyr::summarise(mcnv= mean(mmm_cnv))
+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
@@ -355,15 +357,15 @@ ezANOVA(data = allAverageDat1%>% filter(phase=="Reproduction"), dv= m_cnv, wid=S
 #lmm
 ```{r}
 cnv_build2_r= cnv_build2%>% filter(phase=="Reproduction")
-#cnv_build2_r$targetDur=cnv_build2_r$targetDur/1000
-cnv_build2_r$Dur = cnv_build2_r$targetDur - 1200
+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*Dur+group + # fixed effect, covariate and interaction
-              (Dur|SubName), # random effect
+              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
@@ -382,19 +384,21 @@ mcnv_width_rep %>% group_by(cond) %>% dplyr::summarise(n= n(),time = mean(mlat))
 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$Dur = cnv_width_rep$targetDur - 1200
+cnv_width_rep$targetDur_n = cnv_width_rep$targetDur - 1200
 ```
 
 ```{r}
-scnv_width_rep <- cnv_width_rep %>% group_by(cond, targetDur,Dur, SubName,group) %>% dplyr::summarise(n= n(),time = mean(time))
-#scnv_width_rep$time= scnv_width_rep$time/ 1000
+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*Dur +group+ # fixed effect, covariate and interaction
-              (Dur|SubName), # random effect
+              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
@@ -402,9 +406,9 @@ tab_model(fit_m2_p, p.val = 'kr') # Kenward-Roger approximation, alternative 'sa
 
 # figure 6
 ```{r}
-frow_r <- plot_grid(NULL,p1t_Reproduction_all,nrow=2,ncol = 1, rel_heights  = c(0.35,1), labels=c("b",""))
+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(p_cnv1_target_p,frow_r,nrow=1,ncol = 2, rel_widths  = c(1,4.5), 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

derivatives/behavioral_analysis.Rmd

@@ -1,5 +1,5 @@
 ---
-title: "Exploring range effects in the temporal reproduction using EEG signatures of time"
+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"
 ---
@@ -11,15 +11,22 @@ source('preproc_data.R')
 
 
 ```{r}
-data_beh$rep_err = data_beh$RP - data_beh$targetDur
-
-m_dat= data_beh %>% dplyr::group_by(NSub,targetDur, cond, group) %>%    
+m_dat= data_beh %>% dplyr::group_by(SubName,targetDur, cond, group) %>%    
   dplyr::summarise(n = n(),m_err = mean(rep_err), se_err = sd(rep_err)/ sqrt(n))
 
 mm_dat= m_dat %>% dplyr::group_by(targetDur, cond, group) %>%    
   dplyr::summarise(n = n(),mm_err = mean(m_err), mse_err = sd(m_err)/ sqrt(n))
 ```
 
+```{r}
+m_dat_sub= m_dat %>% dplyr::group_by(SubName, cond, group) %>%    
+  dplyr::summarise(n = n(),mm_err = mean(m_err), mse_err = sd(m_err)/ sqrt(n))
+
+m_dat_sub %>% dplyr::group_by(cond, group) %>%    
+  dplyr::summarise(n = n(),mmm_err = mean(mm_err), mse_err = sd(mm_err)/ sqrt(n))
+```
+
+
 
 ###plot
 ```{r}
@@ -30,8 +37,8 @@ fig_mRep_exp1 = ggplot(mm_dat,
   geom_errorbar(aes(ymin = mm_err- mse_err , ymax = mm_err + mse_err), width = 0.02) + 
   theme_new + scale_color_manual(values = colors_plot) +
   theme(strip.background = element_blank()) +
-  scale_x_continuous(breaks= c(0.4, 0.57, 0.8,1.2,1.7,2.4), labels=c("400","566","800","1200","1700","2400"))+
-  scale_y_continuous(breaks= c(-0.2, -0.1, 0,0.1,0.2), labels=c("-200","-100","0","100","200"))+
+  scale_x_continuous(breaks= c(400, 570, 800,1200,1700,2400), labels=c("400","566","800","1200","1700","2400"))+
+  scale_y_continuous(breaks= c(-200, -100, 0,100,200), labels=c("-200","-100","0","100","200"))+
   labs(x = 'Target interval (ms)', y = 'Reproduction error (ms)', shape= '', color='')+ theme(legend.position=c(0.9,0.95))+ggtitle("")+
   theme(plot.title =element_text(size=12,face='bold')) 
 
@@ -40,23 +47,20 @@ fig_mRep_exp1
 ```
 
 # LMM
-```{r}
-# seconds to miliseconds
-m_dat$targetDur= m_dat$targetDur*1000
-m_dat$m_err= m_dat$m_err*1000
-```
-
 ```{r} 
+m_dat$targetDur_n= m_dat$targetDur-1200 
 m_dat$cond = as.factor(m_dat$cond)
 m_dat$group = as.factor(m_dat$group)
+m_dat= m_dat%>% dplyr::mutate(group= factor(group,levels=c("Upper","Lower")))
 contrasts(m_dat$cond) = contr.Sum(levels(m_dat$cond))
 contrasts(m_dat$group) = contr.Sum(levels(m_dat$group))
-#m_dat$Dur = m_dat$targetDur - 1200
+```
 
+```{r}
 # fit the linear mixed model
-lmm2 = lmer(m_err ~ targetDur+ cond+ group+
-              targetDur:cond + cond:group + # fixed effect and interaction
-              (1|NSub), # random effect
+lmm2 = lmer(m_err ~ targetDur_n+ cond+ group+
+              targetDur_n:cond + cond:group + # fixed effect and interaction
+              (1|SubName), # random effect
               data = m_dat)
 #output as table
 tab_model(lmm2, p.val = 'kr') # Kenward-Roger approximation, alternative 'satterthwaite' approxmiation