library(lme4)
library(lmerTest)
library(confint)
library(emmeans)
library(ggplot2)
library(openxlsx)


dat <- read.xlsx("S3_Appendix_Experimental_data.xlsx",sheet=5)

dat$Participant<-as.factor(dat$Participant)

dat$ExpCond<-as.factor(dat$ExpCond)
dat$Position<-as.factor(dat$Position)
dat$OFB<-as.factor(dat$OFB)
dat$PFB<-as.factor(dat$PFB)

datas<-subset(dat,Position == '-0.5') # use own event onset data after partner's outcome


################## FRN/RewP ##########################
a<-lmer(FRN~ OFB + OFB:( PFB * ExpCond * (SCR + GAR) )+ (1|Participant),data=datas)
summary(a)
anova(a,ddf='Kenward-Roger')
confint(a)


int1=emmeans(a,~OFB|PFB,adjust = "tukey",pbkrtest.limit = 4240,lmerTest.limit = 4240)
pairs(int1,simple="PFB")

p1=emmip(a,OFB~PFB,CIs=TRUE,pbkrtest.limit = 4240,lmerTest.limit = 4240)
p1=p1+labs(y = expression(paste('estimated marginal means [', mu, 'v]')))



int2=emmeans(a,~OFB|ExpCond,adjust = "tukey",pbkrtest.limit = 4240,lmerTest.limit = 4240)
pairs(int2,simple="ExpCond")

p2=emmip(a,OFB~ExpCond,CIs=TRUE,pbkrtest.limit = 4240,lmerTest.limit = 4240)
p2=p2+labs(y = expression(paste('estimated marginal means [', mu, 'v]')))


test(emtrends(a,~ OFB|ExpCond,adjust = "tukey",pbkrtest.limit = 4240,lmerTest.limit = 4240,var ="GAR"))


int3=emtrends(a,~ OFB|ExpCond,var ="GAR",adjust = "tukey",pbkrtest.limit = 4240,lmerTest.limit = 4240)
pairs(int3,simple="ExpCond")



datas2 <- list(GAR=seq(-2,2,by=0.1),OFB=c("correct","incorrect"),ExpCond=c("joint","individual"))
p3=emmip(a,ExpCond~GAR|OFB,at=datas2,CIs=TRUE,pbkrtest.limit = 4240,lmerTest.limit = 4240)
p3=p3+labs(title = 'OFB',y = expression(paste('estimated marginal means [', mu, 'v]')))
p3


p4=ggpubr::ggarrange(p1, p2, p3, labels = LETTERS,nrow=1)
p4



################## P300 ##############################
b<-lmer(P3~ OFB+ OFB : (PFB * ExpCond * (SCR +GAR)) + (1|Participant),data=datas)
summary(b)
anova(b,ddf='Kenward-Roger')
confint(b)


int1=emmeans(b,~OFB|PFB,adjust = "tukey",pbkrtest.limit = 4240,lmerTest.limit = 4240)
pairs(int1,simple="PFB")

p1=emmip(b,OFB~PFB,CIs=TRUE,pbkrtest.limit = 4240,lmerTest.limit = 4240)
p1=p1+labs(y = expression(paste('estimated marginal means [', mu, 'v]')))



int2=emmeans(b,~OFB|ExpCond,adjust = "tukey",pbkrtest.limit = 4240,lmerTest.limit = 4240)
pairs(int2,simple="ExpCond")


p2=emmip(b,OFB~ExpCond,adjust = "tukey",CIs=TRUE,pbkrtest.limit = 4240,lmerTest.limit = 4240)
p2=p2+labs(y = expression(paste('estimated marginal means [', mu, 'v]')))


test(emtrends(b,~ OFB,var ="SCR",adjust = "tukey",pbkrtest.limit = 4240,lmerTest.limit = 4240))


datas2 <- list(SCR=seq(-2,2,by=0.1),OFB=c("correct","incorrect"))
p3=emmip(b,OFB~SCR,at=datas2,CIs=TRUE,pbkrtest.limit = 4240,lmerTest.limit = 4240)
p3=p3+labs(y = expression(paste('estimated marginal means [', mu, 'v]')))


test(emtrends(b,~ OFB|PFB,adjust = "tukey",var ="SCR",pbkrtest.limit = 4240,lmerTest.limit = 4240))


pairs(emtrends(b,~ OFB|PFB,adjust = "tukey",var ="SCR",pbkrtest.limit = 4240,lmerTest.limit = 4240),simple="PFB")


datas2 <- list(SCR=seq(-2,2,by=0.1),OFB=c("correct","incorrect"),PFB=c("correct","incorrect"))
p4=emmip(b,PFB~SCR|OFB,at=datas2,CIs=TRUE,pbkrtest.limit = 4240,lmerTest.limit = 4240)
p4=p4+labs(title = 'OFB',y = expression(paste('estimated marginal means [', mu, 'v]')))

p5=ggpubr::ggarrange(p1, p2, p3, p4, labels = LETTERS)
p5

######################################################