---
title: "CBN reanalysis"
author: "Marijn Struiksma"
date: "08-03-2021"
output: word_document:
    reference_docx: CBN_template.docx
    toc: yes
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, message = FALSE)

```
```{r message=FALSE, warning=FALSE }
library(readr)
library(doBy)
library(psych)
library(lme4)
library(lmerTest)
library(ggplot2)
library(emmeans)
library(sjPlot)
library(sjmisc)
library(sjlabelled)
library(tidyverse)
```

# Analysis procedure

In BrainVision Analyser the data of the CBN study was preprocessed. A Macro copied the trial markers and block markers at the start of the trial (i.e. at the name of the participant) into the segment of the critical word. Baseline correction, artifact rejection on eye-movements (EOG) and on individual channel mode have been performed. The data was also down-sampeled to 500Hz. The trials that contained eye movement artefacts were marked as bad trials. The electrodes that contained artefacts were subsequently marked as bad channels. All critical word segments were exported to mat-files. This procedure was done for 79 participants using the analysis_complete_mat_export.ehtp history file.

A Matlab script (preprocess_trialdata_CIN.m) was used to further process the EEG data in which all trials (segments) of all participants were kept. The script first checked what segments should be removed by looking for the marker 'Bad Interval', gives all segments a name corresponding to the item, and finally replaces the values in bad segments with NA values. Segments are then averaged for the P2 200-250ms (anterior electrodes: Fp1, AF3, AF7, F1, F3, F5, F7, FC1, FC3, FC5, FT7, Fp2, AF4, AF8, F2, F4, F6, F8, FC2,  FC4, FC6, FT8, FPz, AFz, Fz, FCz), N400 350-450ms and LPP 500-650ms (central/posterior electrodes: C1, C3, C5, T7, CP1, CP3, CP5, TP7, P1, P3, P5, P7, PO3, PO7, O1, C2, C4, C6, T8, CP2, CP4, CP6, TP8, P2, P4, P6, P8, PO4, PO8, O2, Cz, CPz, Pz, POz, Oz) to do further trial-based analysis. We've averaged the EEG trial data over max 3 repetitions of the same word per block. 

The SCR data was also exported from BVA in the step mentioned above, and further processed in Matlab (preprocess_trialdata_CIN_GSR.m).
First Ledalab is used to compute the area under the phasic driver. Then the results are selected for the target words (trigger 1-6). We've averaged the SCR trial data over max 3 repetitions of the same word per block. 

The statistical analysis is organized as follows. First I check random part in the multilevel analysis in an empty model. We then start out with a model that contains the three main effects of interest for our hypothesis: Speech Act, Peron and Block (centered). This is our base model. We then stepwise add the main effect of Block squared (centered), followed by 2-way, 3-way and 4-way interactions to end up with the full factorial model. The models are compared using X^2 analysis on the -2LL results. After we've converged to the optimal model we inspect that model and subsequently stepwise remove non-significant main effects.
Finally, we produce a table comparing the optimal model and the full factorial model to answer our hypotheses. The significant effects will be discussed based on the estimates of the optimal final model.

The 4 factors are: Speech Act (Compliment, Neutral, Insult), Person (You, Other), Block Centred (-0.85, 0.85 in 18 steps) and Block2 Centred (the square of Block Centred). For each component (P2, N400 and LPP) we ran separate analyses. 

# Preprocessing

The experiment consisted of 540 trials: 3 Speech Acts, 2 Persons, 3 repetitions of blocks. Each block consisted of 10 unique critical words that were pseudorandomly repeated 3 times within each block. To reduce the noise we've averaged the max 3 repetitions of the same word per block.

Load the EEG data
```{r eval=FALSE}
#setwd("L:/projects/jvbproject/Nicola Project CBN/CBNanalysisTrialdata")
setwd("O:/Research/GW/Projects/Jvbproject/Nicola Project CBN/CBNanalysisTrialdata")
EEG <- read.delim("CBN_ERP_trial_data_13-Nov-2019_79sjs.txt",sep = "\t")
```
```{r eval=FALSE, include = FALSE}
dim(EEG);str(EEG)
```
Average over unique items in each block, max 3 per block, for each participant. Using the mean, excluding NA, and counting the number of items for the average.
```{r eval=FALSE}
sumfun <- function(x, ...){
  c(m=mean(x, na.rm=TRUE, ...), n=length(x))
}
EEGavg<-summaryBy(P2_200.250 + N400_350.450 + LPP_500.650 ~ Participant + block + item, data = EEG,FUN = sumfun)
```
```{r eval=FALSE, include = FALSE}
# checking the averages
dim(EEGavg);str(EEGavg)
head(EEGavg)
hist(EEGavg$P2_200.250.m)
# finding possible outliers in the histogram of the number trials the average is based on
EEGavg[EEGavg$P2_200.250.n>3,]
EEGavg[EEGavg$N400_350.450.n>3,]
EEGavg[EEGavg$LPP_500.650.n>3,]
```

Load the GSR data
```{r eval=FALSE}
GSR <- read.delim("CBN_logGSR_trial_data_04-Jan-2020.txt",sep = ",")
```
```{r eval=FALSE, include = FALSE}
dim(GSR);str(GSR)
```
Average over unique items in each block, max 3 per block, for each participant. Using the mean, excluding NA, and counting the number of items for the average.
```{r eval=FALSE}
GSRavg<-summaryBy(logISCR ~ Participant + block + item, data = GSR,FUN = sumfun)
```
```{r eval=FALSE, include = FALSE}
# checking the averages
dim(GSRavg);str(GSRavg)
head(GSRavg)
hist(GSRavg$logISCR.m,breaks = 30)
# finding possible outliers in the histogram of the number trials the average is based on
GSRavg[GSRavg$logISCR.n>3,]
```
```{r eval=FALSE}
save(EEGavg,GSRavg, file = "CBN_preprocessed_trial_data_avg_79sjs.RData")
```

## Load preprocessed EEG data, averaged per block.
```{r eval=FALSE}
load("CBN_preprocessed_trial_data_avg_79sjs.RData")
EEGdata<-EEGavg
```
```{r eval=FALSE, echo=FALSE, include=FALSE}
#Turn first and third person into a factor variable:
EEGdata$Person<-rep(1,length(EEGdata$item))
EEGdata$Person[EEGdata$item>39]<-3;hist(EEGdata$Person)
EEGdata$Person<-factor(EEGdata$Person,levels = c(1,3),labels = c("1st person","3rd person"))
str(EEGdata$Person)

#Turn the speech act into a factor variable: neutral=2
EEGdata$SpeechAct<-rep(2,length(EEGdata$item))
EEGdata$SpeechAct[EEGdata$item<20]<-1;EEGdata$SpeechAct[39<EEGdata$item&EEGdata$item<50]<-1 #compliment=1
EEGdata$SpeechAct[30<EEGdata$item&EEGdata$item<40]<-3;EEGdata$SpeechAct[59<EEGdata$item]<-3 #insult=3
hist(EEGdata$SpeechAct)
EEGdata$SpeechAct<-factor(EEGdata$SpeechAct,levels = c(1:3),labels = c("compliment","neutral","insult"));str(EEGdata$SpeechAct)

```

## variables:
dependent variables: P2_200.250.m, N400_350.40.m, LPP_500.650.m  
random factors: Participant, item  
independent variables: SpeechAct (1=Compliment, 2=Neutral,3=Insult),Person(1,3),block(continuous, will be centered and divided by 10 for easier interpretation of the estimates)  

```{r eval=FALSE, echo=FALSE, include=FALSE}
#centering block, divide by 10 for easier interpretation of the estimates, and computing block^2
EEGdata$block<-(EEGdata$block-100)/10;head(EEGdata$block)
unique(EEGdata$block)

EEGdata$BlockCent<-EEGdata$block-0.95;unique(EEGdata$BlockCent)
EEGdata$Block2Cent<-EEGdata$BlockCent^2;unique(EEGdata$Block2Cent)
cor(EEGdata$BlockCent,EEGdata$Block2Cent)
```

# 200-250ms window
```{r eval=FALSE, echo=FALSE}
#first check for missing values
which(is.na(EEGdata$P2_200.250.m))
which(is.na(EEGdata$item&EEGdata$P2_200.250.m))
EEGdata<-EEGdata[-6390,]#exclude this row, RUN ONLY ONCE!!!

#check the normal distribution
ggplot(EEGdata,aes(x=P2_200.250.m)) + geom_density() + facet_wrap(Person ~ SpeechAct)

```
```{r eval=FALSE}
save(EEGdata, file = "CBN_EEG_trial_data_for_MLM_analysis_79sjs.RData")
write.table(EEGdata, file = 'CBN_EEG_trial_data_for_MLM_analysis_79sjs.csv',dec = ",",sep = "\t",row.names = FALSE)
```

## Building the models

```{r} 
#load data for analysis, previous steps can be skipped.
load("CBN_EEG_trial_data_for_MLM_analysis_79sjs.RData")
```
```{r eval=TRUE, results="hide"}
#empty model
P2M0<-lmer(P2_200.250.m ~ (1|Participant)+(1|item),data=EEGdata,REML = F);summary(P2M0)
# add 3 main effects
P2M1<-lmer(P2_200.250.m ~ SpeechAct+Person+BlockCent+(1|Participant)+(1|item),data=EEGdata,REML = F);summary(P2M1)
anova(P2M0,P2M1) #p<.001 so P2M1 is a better model
```
```{r eval=TRUE, results="hide"}
# add Block2Cent
P2M2<-lmer(P2_200.250.m ~ SpeechAct+Person+BlockCent+Block2Cent+(1|Participant)+(1|item),data=EEGdata,REML = F);summary(P2M2)
anova(P2M1,P2M2) # p=.07 not significant
# add SpeechAct:Person
P2M3 <- update(P2M2, .~.+SpeechAct:Person);summary(P2M3)
anova(P2M2,P2M3) # p=.75 not significant
# add SpeechAct:BlockCent
P2M4 <- update(P2M3, .~.+SpeechAct:BlockCent);summary(P2M4)
anova(P2M3,P2M4) # p=.88 not significant
# add Person:BlockCent
P2M5 <- update(P2M4, .~.+Person:BlockCent);summary(P2M5)
anova(P2M4,P2M5) # p=.66 not significant
# add SpeechAct:Block2Cent
P2M6 <- update(P2M5, .~.+SpeechAct:Block2Cent);summary(P2M6)
anova(P2M5,P2M6) # p=.86 not significant
# add Person:Block2Cent
P2M7 <- update(P2M6, .~.+Person:Block2Cent);summary(P2M7)
anova(P2M6,P2M7) # p=.70 not significant
# add BlockCent:Block2Cent
P2M8 <- update(P2M7, .~.+BlockCent:Block2Cent);summary(P2M8)
anova(P2M7,P2M8) # p=.26 not significant
# add SpeechAct:Person:Block2Cent
P2M9 <- update(P2M8, .~.+SpeechAct:Person:BlockCent);summary(P2M9)
anova(P2M8,P2M9) # p=.98 not significant
# add SpeechAct:Person:Block2Cent
P2M10 <- update(P2M9, .~.+SpeechAct:Person:Block2Cent);summary(P2M10)
anova(P2M9,P2M10) # p=.33 not significant
# add Person:BlockCent:Block2Cent
P2M11 <- update(P2M10, .~.+Person:BlockCent:Block2Cent);summary(P2M11)
anova(P2M10,P2M11) # p=.07 not significant
# add SpeechAct:Person:BlockCent:Block2Cent
P2M12 <- update(P2M11, .~.+SpeechAct:Person:BlockCent:Block2Cent);summary(P2M12)
anova(P2M11,P2M12) # p=.85 not significant

```
The model P2M1 with the three main effects seems to be the best model. In a next step we compare it to model P2M12.
```{r eval=FALSE}
anova(P2M1,P2M12) # p=.85 so not significant.
```
An alternative model to test is the full-factorial model based on our initial 3 factors, so without Block2Cent.
```{r eval=FALSE}
# full-factorial SpeechAct*Person*BlockCent
P2M1a<-lmer(P2_200.250.m ~ SpeechAct*Person*BlockCent+(1|Participant)+(1|item),data=EEGdata,REML = F);summary(P2M1a)
anova(P2M1,P2M1a) # p=.99 not significant, so P2M1 is a sparser model

```

A final alternative model to test the most sparse model without the non-significant factor Person and only the significant main effects of SpeechAct and BlockCent
```{r eval=TRUE, results="hide"}
# Remove main effect Person
P2M1b<-lmer(P2_200.250.m ~ SpeechAct+BlockCent+(1|Participant)+(1|item),data=EEGdata,REML = F);summary(P2M1b)
anova(P2M1,P2M1b) # p=.19 not significant, so P2M1b is a sparser model

```
```{r eval=FALSE, echo=FALSE}
save(P2M12,P2M1b, file = "CBN_EEG_P2_final_models_79sjs.RData")
load("CBN_EEG_P2_final_models_79sjs.RData")
```
### Summary of best and full factorial Model P2
The best model for the P2 data is P2M1b. We now show the full-factorial model with all 4 factors, but to estimate the parameters correctly the results are based on a model without the non-significant parameters, which is the model with the main effect of SpeechAct and BlockCent.
```{r echo=FALSE, warning=FALSE, P2.lmm.table}
tab_model(P2M1b,P2M12,
          collapse.ci=TRUE, show.stat=TRUE, digits.p=3, ci.hyphen=",&nbsp;", 
          dv.labels=c("P2 best model (M1b)",
                      "** P2 full factorial model (M12) **"),
          show.loglik=TRUE, show.dev=TRUE, show.r2=FALSE, show.obs=FALSE,
          show.df=TRUE, title="Summary of LMM of P2 data." )
tab_model(P2M1b,
          collapse.ci=TRUE, show.stat=TRUE, digits.p=3, ci.hyphen=",&nbsp;", 
          dv.labels=c("P2 best model (M1b)"),
          show.loglik=TRUE, show.dev=TRUE, show.r2=FALSE, show.obs=FALSE,
          show.df=TRUE, title="Summary of LMM of P2 data." )
```
## post hoc tests P2
### Explore SpeechAct
```{r eval=TRUE, warning=FALSE}
#P2M12_emmeans<-emmeans(P2M12,~SpeechAct*Person,pbkrtest.limit = 13903);
P2M1b_emmeans<-emmeans(P2M1b,~SpeechAct);
P2M1b_emmeans
contrast(P2M1b_emmeans,  "tukey",simple = "each", combine = T, adjust = "tukey")
```

### explore SpeechAct:blockCent
```{r eval=TRUE, echo=FALSE}
#P2M1b_emtrends<-emtrends(P2M1b,pairwise~SpeechAct,var="BlockCent",emm_options(pbkrtest.limit = 13903))
P2M1b_emtrends<-emtrends(P2M1b,pairwise~SpeechAct,var="BlockCent")
```
```{r}
P2M1b_emtrends
test(P2M1b_emtrends)
```
```{r eval=TRUE, echo=FALSE, fig.cap = "Figure on SpeechAct 200-250ms with CI"}
#emmip(P2M1b,SpeechAct ~ BlockCent, cov.reduce = range, pbkrtest = 14000, CIs = TRUE)+ scale_color_manual(values=c("#009E73","#E69F00", "#D55E00"))
emmip(P2M1b,SpeechAct ~ BlockCent, cov.reduce = range, CIs = TRUE)+ scale_color_manual(values=c("#009E73","#E69F00", "#D55E00"))
♣```

Plot of regressionlines for the different categories of SpeechAct (slopes don't differ), since there wasn't a significant interaction this figure demonstrates the main effects of SpeechAct and BlockCent:

### Explore SpeechAct:Person
```{r eval=TRUE, warning=FALSE}
#P2M12_emmeans<-emmeans(P2M12,~SpeechAct*Person,pbkrtest.limit = 13903);
P2M12_emmeans<-emmeans(P2M12,~SpeechAct*Person);
P2M12_emmeans
contrast(P2M12_emmeans,  "tukey",simple = "each", combine = T, adjust = "tukey")
```

Since Person is no factor in M1b, this table is generated based on the full factorial model M12. No difference here of course, as Person is no factor in the best model.

Plot of regressionlines for the different categories of SpeechAct:

```{r eval=TRUE, echo=FALSE, fig.cap = "Figure on Person 200-250ms with CI", warning=FALSE}
emmip(P2M12,SpeechAct~Person, CIs = TRUE) + scale_color_manual(values=c("#009E73","#999999", "#D55E00"))
```

## figures P2
```{r, eval=TRUE, echo=FALSE}
#emm_options(pbkrtest.limit = 13903)
P2M1b.df<-emmip(P2M1b,SpeechAct~BlockCent, CIs = TRUE,cov.reduce=F,plotit = F)
```
```{r echo=FALSE}
aa<-ggplot(data=P2M1b.df, aes(x=xvar, y=yvar,group=SpeechAct)) +
  geom_line(linetype = "blank")+
  geom_ribbon(aes(ymin = LCL,ymax = UCL, fill = SpeechAct), alpha = 0.2) + 
  scale_fill_manual(values=c("#009E73","#999999", "#D55E00")) +
  geom_line(aes(colour = SpeechAct), size = 1) +
  scale_color_manual(values=c("#009E73","#999999", "#D55E00")) +
  theme(legend.justification=c(1,1),
        legend.position=c(1,1))+               # Position legend in top right
  #theme(legend.position="left")+
  ggtitle("200-250ms") +
  xlab("BlockCent") + ylab(expression(paste("Mean Voltage (", mu, "V)")))
aa
```


# 350-450ms window
```{r eval=FALSE, echo=FALSE}
#check the normal distribution
ggplot(EEGdata,aes(x=N400_350.450.m)) + geom_density() + facet_wrap(Person ~ SpeechAct)

```
## Building the models

```{r eval=TRUE, results="hide"}
#empty model
N400M0<-lmer(N400_350.450.m ~ (1|Participant)+(1|item),data=EEGdata,REML = F);summary(N400M0)
# add 3 main effects
N400M1<-lmer(N400_350.450.m ~ SpeechAct+Person+BlockCent+(1|Participant)+(1|item),data=EEGdata,REML = F);summary(N400M1)
anova(N400M0,N400M1)
```
```{r eval=TRUE, results="hide"}
# add Block2Cent
N400M2<-lmer(N400_350.450.m ~ SpeechAct+Person+BlockCent+Block2Cent+(1|Participant)+(1|item),data=EEGdata,REML = F);summary(N400M2)
anova(N400M1,N400M2) # p=.03 significant
# add SpeechAct:Person
N400M3 <- update(N400M2, .~.+SpeechAct:Person);summary(N400M3)
anova(N400M2,N400M3) # p=.64 not significant
# add SpeechAct:BlockCent
N400M4 <- update(N400M3, .~.+SpeechAct:BlockCent);summary(N400M4)
anova(N400M3,N400M4) # p=.01 significant
# add Person:BlockCent
N400M5 <- update(N400M4, .~.+Person:BlockCent);summary(N400M5)
anova(N400M4,N400M5) # p=.51 not significant
# add SpeechAct:Block2Cent
N400M6 <- update(N400M5, .~.+SpeechAct:Block2Cent);summary(N400M6)
anova(N400M5,N400M6) # p=.62 not significant
# add Person:Block2Cent
N400M7 <- update(N400M6, .~.+Person:Block2Cent);summary(N400M7)
anova(N400M6,N400M7) # p=.81 not significant
# add BlockCent:Block2Cent
N400M8 <- update(N400M7, .~.+BlockCent:Block2Cent);summary(N400M8)
anova(N400M7,N400M8) # p=.11 not significant
# add SpeechAct:Person:Block2Cent
N400M9 <- update(N400M8, .~.+SpeechAct:Person:BlockCent);summary(N400M9)
anova(N400M8,N400M9) # p=.97 not significant
# add SpeechAct:Person:Block2Cent
N400M10 <- update(N400M9, .~.+SpeechAct:Person:Block2Cent);summary(N400M10)
anova(N400M9,N400M10) # p=.42 not significant
# add Person:BlockCent:Block2Cent
N400M11 <- update(N400M10, .~.+Person:BlockCent:Block2Cent);summary(N400M11)
anova(N400M10,N400M11) # p=.32 not significant
# add SpeechAct:Person:BlockCent:Block2Cent
N400M12 <- update(N400M11, .~.+SpeechAct:Person:BlockCent:Block2Cent);summary(N400M12)
anova(N400M11,N400M12) # p=.97 not significant
```

The model N400M4 with the four main effects seems to be the best model. In a next step we compare it to model N400M12.
```{r eval=FALSE}
anova(N400M4,N400M12) # p=.92 so not significant.
```
Since Person and SpeechAct:Person did not seem to contribute to the model an alternative model to test is a model wihtout the SpeechAct:Person effect.
```{r eval=TRUE, results="hide"}
# Reduced model SpeechAct, Person, BlockCent, Block2Cent, SpeechAct:BlockCent
N400M4a<-lmer(N400_350.450.m ~ SpeechAct + Person + BlockCent + Block2Cent + SpeechAct:BlockCent + (1|Participant)+(1|item),data=EEGdata,REML = F); summary(N400M4a)
anova(N400M4,N400M4a) # p=.65 so not significant, reduced model is sparser

```
```{r eval=FALSE, echo=FALSE}
save(N400M4a,N400M12, file = "CBN_EEG_N400_final_models_79sjs.RData")
load("CBN_EEG_N400_final_models_79sjs.RData")
```
### Summary of best and full factorial Model N400
So removing the SpeechAct:Person interaction results in a sparser model. The best model for the N400 data is N400M4a. We now show the full-factorial model with all 4 factors, but to estimate the parameters correctly the results are based on a model without the non-significant parameters, which is the model with the four main effects and the interaction between SpeechAct and BlockCent.
```{r echo=FALSE, warning=FALSE, N400.lmm.table}
require(sjPlot)
tab_model(N400M4a,N400M12,
          collapse.ci=TRUE, show.stat=TRUE, digits.p=3, ci.hyphen=",&nbsp;", 
          dv.labels=c("N400 best model (M4a)",
                      "** N400 full factorial model (M12) **"),
          show.loglik=TRUE, show.dev=TRUE, show.r2=FALSE, show.obs=FALSE,
          show.df=TRUE, title="Summary of LMM of N400 data." )

tab_model(N400M4a,
          collapse.ci=TRUE, show.stat=TRUE, digits.p=3, ci.hyphen=",&nbsp;", 
          dv.labels=c("N400 best model (M4a)"),
          show.loglik=TRUE, show.dev=TRUE, show.r2=FALSE, show.obs=FALSE,
          show.df=TRUE, title="Summary of LMM of N400 data." )
```
## post hoc tests N400
### Explore SpeechAct
```{r eval=TRUE, warning=FALSE}
#N400M4a_emmeans<-emmeans(N400M4a,~SpeechAct,pbkrtest.limit = 13903);
N400M4a_emmeans<-emmeans(N400M4a,~SpeechAct);
N400M4a_emmeans
contrast(N400M4a_emmeans,  "tukey",simple = "each", combine = T, adjust = "tukey")
```

### explore SpeechAct:blockCent
```{r eval=TRUE, echo=FALSE, warning=FALSE}
#M12aa<-emtrends(M12,pairwise~SpeechAct,var="blockCent",emm_options(pbkrtest.limit = 13903))
N400M4a_emtrends<-emtrends(N400M4a,pairwise~SpeechAct,var="BlockCent")
```
```{r}
N400M4a_emtrends
test(N400M4a_emtrends)
```

The interaction between SpeechAct and BlockCent is significant, the reduction over time differs between compliments and insults. The N400 adapts stronger to insults (-1.31) than to compliments (-0.52).
```{r eval=TRUE, echo=FALSE, fig.cap = "Figure on SpeechAct 350-450ms with CI"}
#emmip(N400M4a,SpeechAct ~ BlockCent, cov.reduce = range, pbkrtest = 14000, CIs = TRUE)+ scale_color_manual(values=c("#009E73","#E69F00", "#D55E00"))
emmip(N400M4a,SpeechAct ~ BlockCent, cov.reduce = range, CIs = TRUE)+ scale_color_manual(values=c("#009E73","#E69F00", "#D55E00"))
```

Plot of regression lines for the different categories of SpeechAct.

### Explore SpeechAct:Person
```{r eval=TRUE, warning=FALSE}
#N400M4a_emmeans<-emmeans(N400M4a,~SpeechAct*Person,pbkrtest.limit = 13903);
N400M4a_emmeans<-emmeans(N400M4a,~SpeechAct*Person);
N400M4a_emmeans
contrast(N400M4a_emmeans,  "tukey",simple = "each", combine = T, adjust = "tukey")
```


Plot of regressionlines for the different categories of SpeechAct:

```{r eval=TRUE, echo=FALSE, fig.cap = "Figure on Person 350-450ms with CI", warning=FALSE}
emmip(N400M4a,SpeechAct~Person, CIs = TRUE)+ scale_color_manual(values=c("#009E73","#999999", "#D55E00"))
```

No difference here of course, the hypotheses as specified compliment first/third person and insult first/third person are not significant in this model (nor in any other). This merely represents the main effects of SpeechAct and Person.

## figures N400
```{r, eval=TRUE, echo=FALSE}
N400M4a.df<-emmip(N400M4a,SpeechAct~BlockCent, CIs = TRUE, cov.reduce=range, plotit = F)
```
```{r echo=FALSE, fig.cap = "Figure on Speech Act over time 350-450ms with CI"}
aa1<-ggplot(data=N400M4a.df, aes(x=xvar, y=yvar,group=SpeechAct)) +
  geom_line(linetype = "blank")+
  geom_ribbon(aes(ymin = LCL,ymax = UCL, fill = SpeechAct), alpha = 0.2) + 
  scale_fill_manual(values=c("#009E73","#999999", "#D55E00")) +
  geom_line(aes(colour = SpeechAct), size = 1) +
  scale_color_manual(values=c("#009E73","#999999", "#D55E00")) +
  theme(legend.justification=c(1,1),
        legend.position=c(1,1))+               # Position legend in top right
  #theme(legend.position="left")+
  ggtitle("N400 350-450ms") +
  xlab("BlockCent") + ylab(expression(paste("Mean Voltage (", mu, "V)")))
aa1
```

# 500-650ms window
```{r eval=FALSE, echo=FALSE}
#check the normal distribution
ggplot(EEGdata,aes(x=LPP_500.650.m)) + geom_density() + facet_wrap(Person ~ SpeechAct)

```
## Building the models

```{r eval=TRUE, results="hide"}
#empty model
LPPM0<-lmer(LPP_500.650.m ~ (1|Participant)+(1|item),data=EEGdata,REML = F);summary(LPPM0)
# add 3 main effects
LPPM1<-lmer(LPP_500.650.m ~ SpeechAct+Person+BlockCent+(1|Participant)+(1|item),data=EEGdata,REML = F);summary(LPPM1)
anova(LPPM0,LPPM1) #p<.001 so M1 is a better model
```
```{r eval=TRUE, results="hide"}
# add Block2Cent
LPPM2<-lmer(LPP_500.650.m ~ SpeechAct+Person+BlockCent+Block2Cent+(1|Participant)+(1|item),data=EEGdata,REML = F);summary(LPPM2)
anova(LPPM1,LPPM2) # p<.001 significant
# add SpeechAct:Person
LPPM3 <- update(LPPM2, .~.+SpeechAct:Person);summary(LPPM3)
anova(LPPM2,LPPM3) # p=.97 not significant
# add SpeechAct:BlockCent
LPPM4 <- update(LPPM3, .~.+SpeechAct:BlockCent);summary(LPPM4)
anova(LPPM3,LPPM4) # p=.001 significant
# add Person:BlockCent
LPPM5 <- update(LPPM4, .~.+Person:BlockCent);summary(LPPM5)
anova(LPPM4,LPPM5) # p=.53 not significant
# add SpeechAct:Block2Cent
LPPM6 <- update(LPPM5, .~.+SpeechAct:Block2Cent);summary(LPPM6)
anova(LPPM5,LPPM6) # p=.96 not significant
# add Person:Block2Cent
LPPM7 <- update(LPPM6, .~.+Person:Block2Cent);summary(LPPM7)
anova(LPPM6,LPPM7) # p=.01  significant
# add BlockCent:Block2Cent
LPPM8 <- update(LPPM7, .~.+BlockCent:Block2Cent);summary(LPPM8)
anova(LPPM7,LPPM8) # p=.64 not significant
# add SpeechAct:Person:Block2Cent
LPPM9 <- update(LPPM8, .~.+SpeechAct:Person:BlockCent);summary(LPPM9)
anova(LPPM8,LPPM9) # p=.44 not significant
# add SpeechAct:Person:Block2Cent
LPPM10 <- update(LPPM9, .~.+SpeechAct:Person:Block2Cent);summary(LPPM10)
anova(LPPM9,LPPM10) # p=.67 not significant
# add Person:BlockCent:Block2Cent
LPPM11 <- update(LPPM10, .~.+Person:BlockCent:Block2Cent);summary(LPPM11)
anova(LPPM10,LPPM11) # p=.04  significant
# add SpeechAct:Person:BlockCent:Block2Cent
LPPM12 <- update(LPPM11, .~.+SpeechAct:Person:BlockCent:Block2Cent);summary(LPPM12)
anova(LPPM11,LPPM12) # p=.45 not significant
```

The model LPPM11 seems to be the best model. In a next step we inspect model LPPM11 and stepwise remove non-significant interactions.
```{r eval=TRUE, results="hide"}
# remove SpeechAct:Person:Block2Cent
LPPM11a<-lmer(LPP_500.650.m ~ SpeechAct + Person + BlockCent + Block2Cent + SpeechAct:Person + SpeechAct:BlockCent + Person:BlockCent + SpeechAct:Block2Cent + Person:Block2Cent + BlockCent:Block2Cent + SpeechAct:Person:BlockCent + Person:BlockCent:Block2Cent + (1 | Participant) + (1 | item), data=EEGdata,REML = F); summary(LPPM11a)
anova(LPPM11a,LPPM11) # p=.68 so not significant, LPPM11a is a sparser model
# remove SpeechAct:Person:BlockCent
LPPM11b<-lmer(LPP_500.650.m ~ SpeechAct + Person + BlockCent + Block2Cent + SpeechAct:Person + SpeechAct:BlockCent + Person:BlockCent + SpeechAct:Block2Cent + Person:Block2Cent + BlockCent:Block2Cent + Person:BlockCent:Block2Cent + (1 | Participant) + (1 | item), data=EEGdata,REML = F); summary(LPPM11b)
anova(LPPM11b,LPPM11a) # p=.42 so not significant, LPPM11b is a sparser model
# remove SpeechAct:Block2Cent
LPPM11c<-lmer(LPP_500.650.m ~ SpeechAct + Person + BlockCent + Block2Cent + SpeechAct:Person + SpeechAct:BlockCent + Person:BlockCent + Person:Block2Cent + BlockCent:Block2Cent + Person:BlockCent:Block2Cent + (1 | Participant) + (1 | item), data=EEGdata,REML = F); summary(LPPM11c)
anova(LPPM11c,LPPM11b) # p=.85 so not significant, LPPM11c is a sparser model
# remove SpeechAct:Person
LPPM11d<-lmer(LPP_500.650.m ~ SpeechAct + Person + BlockCent + Block2Cent + SpeechAct:BlockCent + Person:BlockCent + Person:Block2Cent + BlockCent:Block2Cent + Person:BlockCent:Block2Cent + (1 | Participant) + (1 | item), data=EEGdata,REML = F); summary(LPPM11d)
anova(LPPM11d,LPPM11c) # p=.85 so not significant, LPPM11d is a sparser model
# remove BlockCent:Block2Cent
LPPM11e<-lmer(LPP_500.650.m ~ SpeechAct + Person + BlockCent + Block2Cent + SpeechAct:BlockCent + Person:BlockCent + Person:Block2Cent + Person:BlockCent:Block2Cent + (1 | Participant) + (1 | item), data=EEGdata,REML = F); summary(LPPM11e)
anova(LPPM11e,LPPM11d) # p=1 so not significant, LPPM11e is a sparser model
# remove Person:BlockCent:Block2Cent
LPPM11f<-lmer(LPP_500.650.m ~ SpeechAct + Person + BlockCent + Block2Cent + SpeechAct:BlockCent + Person:BlockCent + Person:Block2Cent + (1 | Participant) + (1 | item), data=EEGdata,REML = F); summary(LPPM11f)
anova(LPPM11f,LPPM11e) # p=.12 so not significant, LPPM11f is a sparser model
# remove Person:BlockCent
LPPM11g<-lmer(LPP_500.650.m ~ SpeechAct + Person + BlockCent + Block2Cent + SpeechAct:BlockCent + Person:Block2Cent + (1 | Participant) + (1 | item), data=EEGdata,REML = F); summary(LPPM11g)
anova(LPPM11g,LPPM11f) # p=.53 so not significant, LPPM11g is a sparser model
```
```{r eval=FALSE, echo=FALSE}
save(LPPM11g,LPPM12, file = "CBN_EEG_LPP_final_models_79sjs.RData")
load("CBN_EEG_LPP_final_models_79sjs.RData")
```
### Summary of best and full factorial Model LPP
Since there are significant interactions the non-significan main effects are kept in the model. The best model for the LPP data is LPPM11g. We now show the full-factorial model with all 4 factors, but to estimate the parameters correctly the results are based on a model without the non-significant parameters, which is the model with the four main effects and the interaction between SpeechAct and BlockCent. The interaction between Person and Block2Cent is also significant.
```{r echo=FALSE, warning=FALSE, LPP.lmm.table}
require(sjPlot)
tab_model(LPPM11g,LPPM12,
          collapse.ci=TRUE, show.stat=TRUE, digits.p=3, ci.hyphen=",&nbsp;", 
          dv.labels=c("LPP best model (M11g)",
                      "** LPP full factorial model (M12) **"),
          show.loglik=TRUE, show.dev=TRUE, show.r2=FALSE, show.obs=FALSE,
          show.df=TRUE, title="Summary of LMM of LPP data." )
tab_model(LPPM11g,
          collapse.ci=TRUE, show.stat=TRUE, digits.p=3, ci.hyphen=",&nbsp;", 
          dv.labels=c("LPP best model (M11g)"),
          show.loglik=TRUE, show.dev=TRUE, show.r2=FALSE, show.obs=FALSE,
          show.df=TRUE, title="Summary of LMM of LPP data." )
```
## post hoc tests LPP
### Explore SpeechAct
```{r eval=TRUE, warning=FALSE}
#LPPM11g_emmeans<-emmeans(LPPM11g,~SpeechAct,pbkrtest.limit = 13903);
LPPM11g_emmeans<-emmeans(LPPM11g,~SpeechAct);
LPPM11g_emmeans
contrast(LPPM11g_emmeans,  "tukey",simple = "each", combine = T, adjust = "tukey")
```

### explore SpeechAct:blockCent
```{r eval=FALSE, warning=FALSE}
load("CBN_EEG_LPP_final_models_79sjs.RData")

```
```{r eval=TRUE, echo=FALSE}
#M103baa<-emtrends(M103b,pairwise~SpeechAct,var="BlockCent",emm_options(pbkrtest.limit = 13903))
LPPM11g_trends<-emtrends(LPPM11g,pairwise~SpeechAct,var="BlockCent")
```
```{r}
LPPM11g_trends
test(LPPM11g_trends)
```
The interaction between SpeechAct and BlockCent is significant, the reduction over time differs between compliments and insults. The LPP adapts stronger to insults (-0.17) than to compliments (-0.07).
```{r eval=TRUE, echo=FALSE, fig.cap = "Figure on SpeechAct 500-650ms with CI's"}
#emmip(LPPM11d,SpeechAct ~ BlockCent, cov.reduce = range, pbkrtest = 14000, CIs = TRUE)+ scale_color_manual(values=c("#009E73","#E69F00", "#D55E00"))
emmip(LPPM11g,SpeechAct ~ BlockCent, cov.reduce = range, CIs = TRUE)+ scale_color_manual(values=c("#009E73","#E69F00", "#D55E00"))
```

Plot of regression lines for the different categories of SpeechAct.

### Explore SpeechAct:Person
```{r eval=TRUE, warning=FALSE}
#LPPM11g+emmeans<-emmeans(LPPM11g,~SpeechAct*Person,pbkrtest.limit = 13903);
LPPM11g_emmeans<-emmeans(LPPM11g,~SpeechAct*Person);
LPPM11g_emmeans
```

```{r}
contrast(LPPM11g_emmeans,  "tukey",simple = "each", combine = T, adjust = "tukey")
```
Plot of regression lines for the different categories of SpeechAct:

```{r eval=TRUE, echo=FALSE, fig.cap = "Figure on Person 500-650ms with CI's", warning=FALSE}
emmip(LPPM11g,SpeechAct~Person, CIs = TRUE)+ scale_color_manual(values=c("#009E73","#999999", "#D55E00"))
```

No difference here of course, the hypotheses as specified compliment first/third person and insult first/third person are not significant in this model (nor in any other).

### explore Person:Block2Cent
```{r}
#M103baa<-emtrends(M103b,pairwise~SpeechAct,var="BlockCent",emm_options(pbkrtest.limit = 13903))
LPPM11g_trends2<-emtrends(LPPM11g,pairwise~Person,var="Block2Cent")
LPPM11g_trends2
```
The interaction between Person and Block2Cent is significant, the reduction over time differs between 1st person and 3rd person. The quadratic component of the LPP is stronger for 1st person (1.69) than 3rd person (0.45).

```{r eval=TRUE, echo=FALSE, fig.cap = "Figure on SpeechAct 500-650ms with CI's"}
#emmip(LPPM11d,SpeechAct ~ BlockCent, cov.reduce = range, pbkrtest = 14000, CIs = TRUE)+ scale_color_manual(values=c("#009E73","#E69F00", "#D55E00"))
emmip(LPPM11g,Person ~ Block2Cent, cov.reduce = range, CIs = TRUE)+ scale_color_manual(values=c("#999999","#000000"))
```

Plot of regression lines for the different categories of Person.

## figures LPP
```{r, eval=TRUE, echo=FALSE}
#emm_options(pbkrtest.limit = 13903)
LPPM11g.df<-emmip(LPPM11g,~SpeechAct~BlockCent, CIs = TRUE,cov.reduce=range,plotit = F)
```
```{r echo=FALSE, fig.cap = "Figure on Speech Act over time 500-650ms"}
aa10<-ggplot(data=LPPM11g.df, aes(x=xvar, y=yvar,group=SpeechAct)) +
  geom_line(linetype = "blank")+
  geom_ribbon(aes(ymin = LCL,ymax = UCL, fill = SpeechAct), alpha = 0.2) + 
  scale_fill_manual(values=c("#009E73","#999999", "#D55E00")) +
  geom_line(aes(colour = SpeechAct), size = 1) +
  scale_color_manual(values=c("#009E73","#999999", "#D55E00")) +
  theme(legend.justification=c(1,1),
        legend.position=c(1,1))+               # Position legend in top right
  #theme(legend.position="left")+
  ggtitle("LPP 500-650ms with CI") +
  xlab("Block centered") + ylab(expression(paste("Mean Voltage (", mu, "V)")))
aa10
```

# SCR analysis
## Load preprocessed SCR data, averaged per block.
```{r eval=FALSE}
load("CBN_preprocessed_trial_data_avg_79sjs.RData")
SCRdata<-GSRavg
```
```{r eval=FALSE, echo=FALSE, include=FALSE}
#Turn first and third person into a factor variable:
SCRdata$Person<-rep(1,length(SCRdata$item))
SCRdata$Person[SCRdata$item>39]<-3;hist(SCRdata$Person)
SCRdata$Person<-factor(SCRdata$Person,levels = c(1,3),labels = c("1st person","3rd person"))
str(SCRdata$Person)

#Turn the speech act into a factor variable: neutral=2
SCRdata$SpeechAct<-rep(2,length(SCRdata$item))
SCRdata$SpeechAct[SCRdata$item<20]<-1;SCRdata$SpeechAct[39<SCRdata$item&SCRdata$item<50]<-1 #compliment=1
SCRdata$SpeechAct[30<SCRdata$item&SCRdata$item<40]<-3;SCRdata$SpeechAct[59<SCRdata$item]<-3 #insult=3
hist(SCRdata$SpeechAct)
SCRdata$SpeechAct<-factor(SCRdata$SpeechAct,levels = c(1:3),labels = c("compliment","neutral","insult"));str(SCRdata$SpeechAct)

```

## variables:
dependent variable: logISCR.m  
random factors: Participant, item  
independent variables: SpeechAct (1=Compliment, 2=Neutral,3=Insult),Person(1,3),block(continuous, will be centered and divided by 10 for easier interpretation of the estimates)  

```{r eval=FALSE, echo=FALSE, include=FALSE}
#centering block, divide by 10 for easier interpretation of the estimates, and computing block^2
SCRdata$block<-(SCRdata$block-100)/10;head(SCRdata$block)
unique(SCRdata$block)

SCRdata$BlockCent<-SCRdata$block-0.95;unique(SCRdata$BlockCent)
SCRdata$Block2Cent<-SCRdata$BlockCent^2;unique(SCRdata$Block2Cent)
cor(SCRdata$BlockCent,SCRdata$Block2Cent)
```
```{r eval=FALSE, echo=FALSE}
#first check for missing values
which(is.na(SCRdata$logISCR.m))
which(is.na(SCRdata$item&SCRdata$logISCR.m))

#check the normal distribution
ggplot(SCRdata,aes(x=logISCR.m)) + geom_density() + facet_wrap(Person ~ SpeechAct)

```

```{r eval=FALSE}
save(SCRdata, file = "CBN_SCR_trial_data_for_MLM_analysis_79sjs.RData")
write.table(SCRdata, file = 'CBN_SCR_trial_data_for_MLM_analysis_79sjs.csv',dec = ",",sep = "\t",row.names = FALSE)
```

## Building the models

```{r} 
#load data for analysis, previous steps can be skipped.
load("CBN_SCR_trial_data_for_MLM_analysis_79sjs.RData")
```
```{r eval=TRUE, results="hide"}
#empty model
SCRM0<-lmer(logISCR.m ~ (1|Participant)+(1|item),data=SCRdata,REML = F);summary(SCRM0)
# add 3 main effects
SCRM1<-lmer(logISCR.m ~ SpeechAct+Person+BlockCent+(1|Participant)+(1|item),data=SCRdata,REML = F);summary(SCRM1)
# SCRM1 is singular, the item variance is close to 0. An alternative is to remove item from the random part
SCRM1a<-lmer(logISCR.m ~ SpeechAct+Person+BlockCent+(1|Participant),data=SCRdata,REML = F);summary(SCRM1a)
anova(SCRM0,SCRM1a) #p<.001 so P2M1a is a better model
```

```{r eval=TRUE, results="hide"}
# add Block2Cent
SCRM2<-lmer(logISCR.m~ SpeechAct+Person+BlockCent+Block2Cent+(1|Participant),data=SCRdata,REML = F);summary(SCRM2)
anova(SCRM1a,SCRM2) # p<.001 significant
# add SpeechAct:Person
SCRM3 <- update(SCRM2, .~.+SpeechAct:Person);summary(SCRM3)
anova(SCRM2,SCRM3) # p=.63 not significant
# add SpeechAct:BlockCent
SCRM4 <- update(SCRM3, .~.+SpeechAct:BlockCent);summary(SCRM4)
anova(SCRM3,SCRM4) # p=.70 not significant
# add Person:BlockCent
SCRM5 <- update(SCRM4, .~.+Person:BlockCent);summary(SCRM5)
anova(SCRM4,SCRM5) # p=.001 significant
# add SpeechAct:Block2Cent
SCRM6 <- update(SCRM5, .~.+SpeechAct:Block2Cent);summary(SCRM6)
anova(SCRM5,SCRM6) # p=.14 not significant
# add Person:Block2Cent
SCRM7 <- update(SCRM6, .~.+Person:Block2Cent);summary(SCRM7)
anova(SCRM6,SCRM7) # p<.001 significant
# add BlockCent:Block2Cent
SCRM8 <- update(SCRM7, .~.+BlockCent:Block2Cent);summary(SCRM8)
anova(SCRM7,SCRM8) # p<.001 significant
# add SpeechAct:Person:Block2Cent
SCRM9 <- update(SCRM8, .~.+SpeechAct:Person:BlockCent);summary(SCRM9)
anova(SCRM8,SCRM9) # p=.004 significant
# add SpeechAct:Person:Block2Cent
SCRM10 <- update(SCRM9, .~.+SpeechAct:Person:Block2Cent);summary(SCRM10)
anova(SCRM9,SCRM10) # p=.14 not significant
# add Person:BlockCent:Block2Cent
SCRM11 <- update(SCRM10, .~.+Person:BlockCent:Block2Cent);summary(SCRM11)
anova(SCRM10,SCRM11) # p=.052 not significant
# add SpeechAct:Person:BlockCent:Block2Cent
SCRM12 <- update(SCRM11, .~.+SpeechAct:Person:BlockCent:Block2Cent);summary(SCRM12)
anova(SCRM11,SCRM12) # p=.056 not significant
```

The model SCRM9 seems to be the best model. In a next step we inspect model SCRM9 and stepwise remove non-significant interactions.
```{r eval=TRUE, results="hide"}
# remove SpeechAct:Person:Block2Cent
SCRM9a<-lmer(logISCR.m ~ SpeechAct + Person + BlockCent + Block2Cent + SpeechAct:Person + SpeechAct:BlockCent + Person:BlockCent + SpeechAct:Block2Cent + Person:Block2Cent + BlockCent:Block2Cent + (1 | Participant), data=SCRdata,REML = F); summary(SCRM9a)
anova(SCRM9a,SCRM9) # p=.68 so not significant, SCRM9a is a sparser model
# remove SpeechAct:Block2Cent
SCRM9b<-lmer(logISCR.m ~ SpeechAct + Person + BlockCent + Block2Cent + SpeechAct:Person + SpeechAct:BlockCent + Person:BlockCent + Person:Block2Cent + BlockCent:Block2Cent + (1 | Participant), data=SCRdata,REML = F); summary(SCRM9b)
anova(SCRM9b,SCRM9a) # p=.20 so not significant, SCRM9b is a sparser model
# remove SpeechAct:BlockCent
SCRM9c<-lmer(logISCR.m ~ SpeechAct + Person + BlockCent + Block2Cent + SpeechAct:Person + Person:BlockCent + Person:Block2Cent + BlockCent:Block2Cent + (1 | Participant), data=SCRdata,REML = F); summary(SCRM9c)
anova(SCRM9c,SCRM9b) # p=.73 so not significant, SCRM9c is a sparser model
# remove SpeechAct:Person
SCRM9d<-lmer(logISCR.m ~ SpeechAct + Person + BlockCent + Block2Cent + Person:BlockCent + Person:Block2Cent + BlockCent:Block2Cent + (1 | Participant), data=SCRdata,REML = F); summary(SCRM9d)
anova(SCRM9d,SCRM9c) # p=.62 so not significant, SCRM9d is a sparser model
```
```{r eval=FALSE, echo=FALSE}
save(SCRM9d,SCRM12, file = "CBN_EEG_SCR_final_models_79sjs.RData")
load("CBN_EEG_SCR_final_models_79sjs.RData")
```
### Summary of best and full factiorial Model SCR
Since there are significant interactions the non-significan main effects are kept in the model. The best model for the SCR data is SCRM9d. We now show the full-factorial model with all 4 factors, but to estimate the parameters correctly the results are based on a model without the non-significant parameters, which is the model with the four main effects and the interaction between Person and BlockCent, the interaction between Person and Block2Cent and the interaction between BlockCent and Block2Cent.
```{r echo=FALSE, warning=FALSE, SCR.lmm.table}
require(sjPlot)
tab_model(SCRM9d,SCRM12,
          collapse.ci=TRUE, show.stat=TRUE, digits.p=3, ci.hyphen=",&nbsp;", 
          dv.labels=c("SCR best model (M9d)",
                      "** SCR full factorial model (M12) **"),
          show.loglik=TRUE, show.dev=TRUE, show.r2=FALSE, show.obs=FALSE,
          show.df=TRUE, title="Summary of LMM of SCR data." )
tab_model(SCRM9d,
          collapse.ci=TRUE, show.stat=TRUE, digits.p=3, ci.hyphen=",&nbsp;", 
          dv.labels=c("SCR best model (M9d)"),
          show.loglik=TRUE, show.dev=TRUE, show.r2=FALSE, show.obs=FALSE,
          show.df=TRUE, title="Summary of LMM of SCR data." )
```
## post hoc tests SCR
### Explore SpeechAct
```{r eval=TRUE, warning=FALSE}
#SCRM9d_emmeans<-emmeans(SCRM9d,~SpeechAct,pbkrtest.limit = 13903);
SCRM9d_emmeans<-emmeans(SCRM9d,~SpeechAct);
SCRM9d_emmeans
contrast(SCRM9d_emmeans,  "tukey",simple = "each", combine = T, adjust = "tukey")
```

### explore Person:BlockCent
```{r eval=FALSE, warning=FALSE}
load("CBN_EEG_SCR_final_models_79sjs.RData")

```
```{r}
SCRM9d_trends<-emtrends(SCRM9d,pairwise~Person,var="BlockCent")
SCRM9d_trends
```
The interaction between Person and BlockCent is significant, the reduction over time differs between 1st and 3rd person. The SCR adapts stronger to 1st person (-0.043) than to compliments (-0.025).
```{r eval=TRUE, echo=FALSE, fig.cap = "Figure on Person:BlockCent SCR with CI's"}
emmip(SCRM9d,Person ~ BlockCent, cov.reduce = range, CIs = TRUE)+ scale_color_manual(values=c("#999999","#000000"))
```

Plot of regression lines for the different categories of Person.

```{r, eval=TRUE, echo=FALSE}
#emm_options(pbkrtest.limit = 13903)
SCR9Md.df<-emmip(SCRM9d,~Person~BlockCent, CIs = TRUE,cov.reduce=range,plotit = F)
```
```{r echo=FALSE, fig.cap = "Figure on Person over time for SCR"}
aa20<-ggplot(data=SCR9Md.df, aes(x=xvar, y=yvar,group=Person)) +
  geom_line(linetype = "blank")+
  geom_ribbon(aes(ymin = LCL,ymax = UCL, fill = Person), alpha = 0.2) + 
  scale_fill_manual(values=c("#999999","#000000")) +
  geom_line(aes(colour = Person), size = 1) +
  scale_color_manual(values=c("#999999","#000000")) +
  theme(legend.justification=c(1,1),
        legend.position=c(1,1))+               # Position legend in top right
  #theme(legend.position="left")+
  ggtitle("SCR with CI") +
  xlab("Block centered") + ylab(expression(paste("Mean Voltage (", mu, "V)")))
aa20
```

### explore Person:Block2Cent
```{r}
SCRM9d_trends2<-emtrends(SCRM9d,pairwise~Person,var="Block2Cent")
SCRM9d_trends2
```
The interaction between Person and Block2Cent is significant, the quadratic reduction over time differs between 1st and 3rd person. The SCR adapts stronger to 1st person (0.101) than to compliments (0.037).
```{r eval=TRUE, echo=FALSE, fig.cap = "Figure on Person:Block2Cent SCR with CI's"}
emmip(SCRM9d,Person ~ Block2Cent, cov.reduce = range, CIs = TRUE)+ scale_color_manual(values=c("#999999","#000000"))
```

Plot of regression lines for the different categories of Person.

### explore SpeechAct:Person
```{r eval=TRUE, warning=FALSE}
#P2M12_emmeans<-emmeans(P2M12,~SpeechAct*Person,pbkrtest.limit = 13903);
SCRM9d_emmeans<-emmeans(SCRM9d,~SpeechAct*Person);
SCRM9d_emmeans
contrast(SCRM9d_emmeans,  "tukey",simple = "each", combine = T, adjust = "tukey")
```

Plot of regressionlines for the different categories of SpeechAct:

```{r eval=TRUE, echo=FALSE, fig.cap = "Figure on Person 200-250ms with CI", warning=FALSE}
emmip(SCRM9d,SpeechAct~Person, CIs = TRUE)+ scale_color_manual(values=c("#009E73","#999999", "#D55E00"))
```

## figures SpeechAct*BlockCent

```{r, eval=TRUE, echo=FALSE}
#emm_options(pbkrtest.limit = 13903)
SCRM9d<-emmip(SCRM9d,~SpeechAct~BlockCent, CIs = TRUE,cov.reduce=range,plotit = F)
```
```{r echo=FALSE, fig.cap = "Figure on Speech Act over time 500-650ms"}
aa20<-ggplot(data=SCRM9d, aes(x=xvar, y=yvar,group=SpeechAct)) +
  geom_line(linetype = "blank")+
  geom_ribbon(aes(ymin = LCL,ymax = UCL, fill = SpeechAct), alpha = 0.2) + 
  scale_fill_manual(values=c("#009E73","#999999", "#D55E00")) +
  geom_line(aes(colour = SpeechAct), size = 1) +
  scale_color_manual(values=c("#009E73","#999999", "#D55E00")) +
  theme(legend.justification=c(1,1),
        legend.position=c(1,1))+               # Position legend in top right
  #theme(legend.position="left")+
  ggtitle("SCR with CI") +
  xlab("Block centered") + ylab(expression(paste("log ISCR (", mu, "S^2)")))
aa20
```

## figures predicted values

In order to interpret the three 2-way interactions we creat a simplified factor of Block*Block2 and plot the results of Person against the development over time.

```{r, eval=TRUE, echo=FALSE}
library(tidyverse)
library(ggeffects)
```
``` {r, eval=TRUE, echo=FALSE}
# create simplified Condition factor for plotting
SCRdata_plot_1st <- SCRdata %>% 
  filter(Person == "1st person")%>% 
  mutate(Condition = case_when(SpeechAct=="compliment" ~ "1st person compliment",
                        SpeechAct=="neutral" ~ "1st person neutral",
                        SpeechAct=="insult" ~ "1st person insult"))
SCRdata_plot_3rd <- SCRdata %>% 
  filter(Person == "3rd person")%>% 
  mutate(Condition = case_when(SpeechAct=="compliment" ~ "3rd person compliment",
                        SpeechAct=="neutral" ~ "3rd person neutral",
                        SpeechAct=="insult" ~ "3rd person insult"))

SCRdata_plot <- bind_rows(SCRdata_plot_1st, SCRdata_plot_3rd) 
SCRdata_plot <- SCRdata_plot %>% 
  mutate(time = BlockCent+Block2Cent)
SCRdata_plot$Condition <- as.factor(SCRdata_plot$Condition)
```
```{r, eval=TRUE, echo=FALSE}
SCRdata_plot$BlockBlock2<-SCRdata_plot$BlockCent*SCRdata_plot$Block2Cent 
SCRM9d_plot4<-lmer(logISCR.m ~ SpeechAct + BlockBlock2*Person
                + (1 | Participant), 
             data=SCRdata_plot,REML = F); #summary(SCRM9d_plot4)
require(sjPlot)
plot_model(SCRM9d_plot4, type = "int",  ci.lvl = 0.95,mdrt.values = "zeromax")
```

```{r, eval=FALSE, echo=FALSE}
SCRM9d_plot<-lmer(logISCR.m ~ Condition + Person + BlockCent + Block2Cent + Person:BlockCent + Person:Block2Cent + BlockCent:Block2Cent + (1 | Participant), data=SCRdata_plot,REML = F)

pred_SCRM9d_plot <- ggpredict(SCRM9d_plot, terms = c("Block2Cent", "BlockCent", "Condition", "Person"), type = "fixed")
plot(pred_SCRM9d_plot)

SCRM9d_plot2<-lmer(logISCR.m ~ Condition + time + Person:time + (1 | Participant), data=SCRdata_plot,REML = F)

pred_SCRM9d_plot2 <- ggpredict(SCRM9d_plot2, terms = c("time", "Condition"), type = "fixed")
plot(pred_SCRM9d_plot2)

SCRM9d_plot3<-lmer(logISCR.m ~ SpeechAct + Person + time + Person:time + (1 | Participant), data=SCRdata_plot,REML = F)

pred_SCRM9d_plot3 <- ggpredict(SCRM9d_plot3, terms = c("time", "Person"), type = "fixed")
plot(pred_SCRM9d_plot3)

SCRM9d_plot5<-lmer(logISCR.m ~ SpeechAct + Person + BlockCent + Block2Cent + 
      Person:BlockCent + Person:Block2Cent + BlockCent:Block2Cent + (1 | Participant), 
      data=SCRdata,REML = F); summary(SCRM9d)
plot_model(SCRM9d, type = "int",  ci.lvl = 0.95,mdrt.values = "zeromax")

```

## Figures main effect of Speech Act (data)

```{r, eval=FALSE, echo=FALSE}
# create 3 blocks for plotting and compute 95% CI
SCRdata_plot_SA <- SCRdata %>% 
  group_by(SpeechAct)  %>% 
  summarise(SCRmean = mean(logISCR.m)[1],
            SCRci = ci(logISCR.m)[1])

colorset = c("#009E73","#999999", "#D55E00")
aa21 <- ggplot(SCRdata_plot_SA) +
  geom_col( aes(x=SpeechAct, y=SCRmean), fill = colorset, alpha=0.5) +
  ylim(0,0.31) +
  geom_errorbar( aes(x=SpeechAct, ymin=SCRmean-SCRci, ymax=SCRmean+SCRci), colour = colorset, width=0.4, alpha=0.9, size=1.0)  +
  xlab("Speech Act") + ylab(expression(paste("log ISCR (", mu, "S^2)"))) +
  coord_flip()
aa21
```

## Figures interaction of Speech Act per block (data)

```{r, eval=FALSE, echo=FALSE}
# create 3 blocks for plotting and compute 95% CI
SCRdata_plot_SAB <- SCRdata %>% 
  mutate(Block3 = case_when(block < 0.65 ~ 1,
                        block < 1.25 ~ 2,
                        block > 1.25 ~ 3))%>% 
  mutate(SpeechAct_num = case_when(SpeechAct=="compliment" ~ 10,
                                        SpeechAct=="neutral" ~ 20,
                                        SpeechAct=="insult"~30)) %>% 
  mutate(SpeechAct_by_Block = SpeechAct_num+Block3) %>% 
  group_by(SpeechAct_by_Block)  %>% 
  summarise(Block3 = first(Block3),
            SpeechAct = first(SpeechAct),
            SCRmean = mean(logISCR.m)[1],
            SCRci = ci(logISCR.m)[1])

# Plot results from block 1
SCRdata_plot_SAB_1 <- SCRdata_plot_SAB %>% 
  filter(Block3 == 1)
  
colorset = c("#009E73","#999999", "#D55E00")
aa22 <- ggplot(SCRdata_plot_SAB_1) +
  geom_col( aes(x=SpeechAct, y=SCRmean), fill = colorset, alpha=0.5) +
  ylim(0,0.31) +
  geom_errorbar( aes(x=SpeechAct, ymin=SCRmean-SCRci, ymax=SCRmean+SCRci), colour = colorset, width=0.4, alpha=0.9, size=1.0)  +
  xlab("Speech Act") + ylab(expression(paste("log ISCR (", mu, "S^2)"))) +
  coord_flip()
aa22

# Plot results from block 2
SCRdata_plot_SAB_2 <- SCRdata_plot_SAB %>% 
  filter(Block3 == 2)
  
aa23 <- ggplot(SCRdata_plot_SAB_2) +
  geom_col( aes(x=SpeechAct, y=SCRmean), fill = colorset, alpha=0.5) +
  ylim(0,0.31) +
  geom_errorbar( aes(x=SpeechAct, ymin=SCRmean-SCRci, ymax=SCRmean+SCRci), colour = colorset, width=0.4, alpha=0.9, size=1.0)  +
  xlab("Speech Act") + ylab(expression(paste("log ISCR (", mu, "S^2)"))) +
  coord_flip()
aa23

# Plot results from block 3
SCRdata_plot_SAB_3 <- SCRdata_plot_SAB %>% 
  filter(Block3 == 3)
  
aa24 <- ggplot(SCRdata_plot_SAB_3) +
  geom_col( aes(x=SpeechAct, y=SCRmean), fill = colorset, alpha=0.5) +
  ylim(0,0.31) +
  geom_errorbar( aes(x=SpeechAct, ymin=SCRmean-SCRci, ymax=SCRmean+SCRci), colour = colorset, width=0.4, alpha=0.9, size=1.0)  +
  xlab("Speech Act") + ylab(expression(paste("log ISCR (", mu, "S^2)"))) +
  coord_flip()
aa24
```

## Figures interaction of Speech Act and Person (data)

```{r, eval=FALSE, echo=FALSE}
# create 3 blocks for plotting and compute 95% CI
z_star_95 <- qnorm(0.975)

SCRdata_plot_SAP <- SCRdata %>% 
  mutate(Person_num = case_when(Person == "1st person" ~ 1,
                        Person == "3rd person" ~ 3)) %>% 
  mutate(SpeechAct_num = case_when(SpeechAct=="compliment" ~ 10,
                                        SpeechAct=="neutral" ~ 20,
                                        SpeechAct=="insult"~30)) %>% 
  mutate(SpeechAct_by_Person = SpeechAct_num+Person_num) %>% 
  group_by(SpeechAct_by_Person)  %>% 
  summarise(Person = first(Person),
            SpeechAct = first(SpeechAct),
            SCRmean = mean(logISCR.m),
            SCRnum = n_distinct(logISCR.m),
            SCRlower = mean(logISCR.m) - z_star_95 * (sd(logISCR.m) / sqrt(SCRnum)),
            SCRupper = mean(logISCR.m) + z_star_95 * (sd(logISCR.m) / sqrt(SCRnum)))
            
# Plot results from 1st Person You
SCRdata_plot_SAP_1st <- SCRdata_plot_SAP %>% 
  filter(Person == "1st person")
  
colorset = c("#009E73","#999999", "#D55E00")
aa25 <- ggplot(SCRdata_plot_SAP_1st) +
  geom_col( aes(x=SpeechAct, y=SCRmean), fill = colorset, alpha=0.5) +
  ylim(0,0.31) +
  geom_errorbar( aes(x=SpeechAct, ymin=SCRlower, ymax=SCRupper), colour = colorset, width=0.4, alpha=0.9, size=1.0)  +
  xlab("Speech Act") + ylab(expression(paste("log ISCR (", mu, "S^2)"))) +
  coord_flip()
aa25

# Plot results from 1st Person You
SCRdata_plot_SAP_3rd <- SCRdata_plot_SAP %>% 
  filter(Person == "3rd person")

aa26 <- ggplot(SCRdata_plot_SAP_3rd) +
  geom_col( aes(x=SpeechAct, y=SCRmean), fill = colorset, alpha=0.5) +
  ylim(0,0.31) +
  geom_errorbar( aes(x=SpeechAct, ymin=SCRlower, ymax=SCRupper), colour = colorset, width=0.4, alpha=0.9, size=1.0)  +
  xlab("Speech Act") + ylab(expression(paste("log ISCR (", mu, "S^2)"))) +
  coord_flip()
aa26
```