trial_level_analysis.Rmd

---
title: "Trial-level analysis"
output: html_document
date: "2023-10-18"
editor_options: 
  chunk_output_type: console
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(plyr)
library(mgcv)
library(itsadug)
library(RColorBrewer)
library(R.matlab)
library(vwr)
library(corrplot)
# Color for conditions
words_col <- brewer.pal(12, 'Paired')[6]
non_word_col <- c('#17bce5',brewer.pal(12, 'Paired')[2])

# Color for recovered stages
colrc <- brewer.pal(7, "Paired")
```

## Data Loading
```{r}
# Collect Trial-data and stimulus info for every trial
hsmm_dat <- read.csv("./results/trial_dat_hsmm_final_excl_6_15_bs_100_.csv")
stimuli <- read.csv("./data/stimuli.csv",sep=";")

stimuli$logSubtlexW <- NA
stimuli$logCelexW <- NA
stimuli$logOLDW <- NA
stimuli$logOLDPW <- NA
stimuli$logOLDRS <- NA

# Load in dlp item-level data
load("./data/dlp-items.Rdata")
load("./data/dlp-stimuli.Rdata")
dlp <- merge(dlp.items,dlp.stimuli)
rm(dlp.items,dlp.stimuli)

# Extract words for easy OLD calculation
dlp_words <- as.character(droplevels(dlp$spelling[dlp$lexicality == "W"]))

# Subtlex vs. Celex differences example:
sum((dlp$subtlex.frequency.million[!is.na(dlp$subtlex.frequency.million)] > dlp$subtlex.frequency.million[dlp$spelling == "vermoord"])*1)/nrow(dlp[!is.na(dlp$subtlex.frequency.million),])
sum((dlp$celex.frequency[!is.na(dlp$celex.frequency)] > dlp$celex.frequency[dlp$spelling == "vermoord"])*1)/nrow(dlp[!is.na(dlp$celex.frequency),])

hsmm_dat$subjects_var <- as.factor(hsmm_dat$subjects_var)
hsmm_dat$blocks_var <- as.factor(hsmm_dat$blocks_var)

hsmm_dat$conditions[hsmm_dat$conditions == 1] <- "W"
hsmm_dat$conditions[hsmm_dat$conditions == 2] <- "PW"
hsmm_dat$conditions[hsmm_dat$conditions == 3] <- "RS"
hsmm_dat$conditions <- as.factor(hsmm_dat$conditions)

hsmm_dat$logGoogleFreq <- NA
hsmm_dat$foreign <- NA
hsmm_dat$logSubtlex <- NA
hsmm_dat$logOLD <- NA

for(stim in unique(hsmm_dat$stims_var)){
  stim_cond <- unique(hsmm_dat$conditions[hsmm_dat$stims_var == stim])
  
  if(stim_cond == "W"){
    hsmm_dat$logGoogleFreq[hsmm_dat$stims_var == stim] <- stimuli$logFreqW[stimuli$word == stim]
    hsmm_dat$foreign[hsmm_dat$stims_var == stim] <- stimuli$foreign[stimuli$word == stim]
    
    # Subtlex per million, taken from DLP corpus (Keuleers et al., 2010) and log-transformed
    hsmm_dat$logSubtlex[hsmm_dat$stims_var == stim] <- log(dlp$subtlex.frequency.million[dlp$spelling == stim])
    stimuli$logSubtlexW[stimuli$word == stim] <- log(dlp$subtlex.frequency.million[dlp$spelling == stim])
    
    # log(OLD20) for stimulus against all words in the DLP (see Yarkoni et al., 2008; Hendrix & Sun, 2020).
    stimuli$logOLDW[stimuli$word == stim] <- log(old20(stim,dlp_words))
    
    # Celex frequency (Baayen et al., 1993) also taken from DLP.
    stimuli$logCelexW[stimuli$word == stim] <- log(dlp$celex.frequency[dlp$spelling == stim])
  }
  
  else if(stim_cond == "PW"){
    hsmm_dat$logGoogleFreq[hsmm_dat$stims_var == stim] <- stimuli$logFreqPW[stimuli$pseudo == stim]
    stimuli$logOLDPW[stimuli$pseudo == stim] <- log(old20(stim,dlp_words))
  }
  
  else if(stim_cond == "RS"){
    hsmm_dat$logGoogleFreq[hsmm_dat$stims_var == stim] <- stimuli$logFreqRS[stimuli$random == stim]
    stimuli$logOLDRS[stimuli$random == stim] <- log(old20(stim,dlp_words))
  }
  
  else {
    stop("Stim not found")
  }
}
stimuli$logCelexW[!is.finite(stimuli$logCelexW)] <- NA

# Collect trial-level EEG data as well
matlab_dat <- readMat("./data/data4HMM_baseline_LD_excl_6_15_bs_100_.mat")

# Combine shared fifth stage and extra stage for pseudo-words
hsmm_dat$durations_65 <- hsmm_dat$durations_6 + hsmm_dat$durations_5
```

## Correlation matrix Google frequency
```{r}
w_pred_dat <- stimuli[stimuli$foreign == 0,colnames(stimuli) %in% c("logFreqW","logSubtlexW","logOLDW","logCelexW")]
lw_pred_dat <- stimuli[stimuli$foreign == 1,colnames(stimuli) %in% c("logFreqW","logSubtlexW","logOLDW","logCelexW")]
pw_pred_dat <- stimuli[stimuli$pseudo != "",colnames(stimuli) %in% c("logFreqPW","logSubtlexW","logOLDPW","logCelexW")]
rs_pred_dat <- stimuli[stimuli$random != "",colnames(stimuli) %in% c("logFreqRS","logSubtlexW","logOLDRS","logCelexW")]

colnames(w_pred_dat) <- c("Google","Subtlex","Celex","OLD20")
colnames(lw_pred_dat) <- c("Google","Subtlex","Celex","OLD20")
colnames(pw_pred_dat) <- c("Google","Subtlex","Celex","OLD20")
colnames(rs_pred_dat) <- c("Google","Subtlex","Celex","OLD20")

cor(w_pred_dat,w_pred_dat,use = "pairwise")
cor(lw_pred_dat,lw_pred_dat,use = "pairwise")
cor(pw_pred_dat,pw_pred_dat,use = "pairwise")
cor(rs_pred_dat,rs_pred_dat,use = "pairwise")

```

## Google Frequency Overview
```{r}
par(mfrow=c(1,1))

# Correlation plots
plot(stimuli$logFreqW,stimuli$logSubtlexW,
     col=words_col,pch=16,
     xlab="Log(Google frequency)",
     ylab="Log(SUBTLEX frequency)",
     xlim=c(0,25))

# Separate color for foreign words
points(stimuli$logFreqW[stimuli$foreign == 1],
       stimuli$logSubtlexW[stimuli$foreign == 1],
       col="purple",pch=16)

plot(stimuli$logFreqPW,stimuli$logSubtlexW,
     col=non_word_col[2],pch=16,
     xlab="Log(Google frequency)",
     ylab="Log(SUBTLEX frequency)",
     xlim=c(0,25))

plot(stimuli$logFreqRS,stimuli$logOLDRS,
     col=non_word_col[1],pch=16,
     xlab="Log(Google frequency)",
     ylab="Log(OLD20)",
     xlim=c(0,25))

# Densities
emptyPlot(xlim=c(0,25),ylim=c(0,0.2),
       xlab="Log(Google frequency)",
       ylab="Density",main="")

lines(density(stimuli$logFreqW[stimuli$foreign == 1],cut=0),col = alpha("purple",0.6),lwd=2)
lines(density(stimuli$logFreqW[stimuli$foreign == 0],cut=0),col = alpha(words_col,1),lwd=2)
lines(density(stimuli$logFreqPW[!is.na(stimuli$logFreqPW)],cut=0),col = alpha(non_word_col[2],1),lwd=2)
lines(density(stimuli$logFreqRS[!is.na(stimuli$logFreqRS)],cut=0),col = alpha(non_word_col[1],1),lwd=2)

emptyPlot(xlim=c(0,1.5),ylim=c(0,6),
       xlab="Log(OLD20)",
       ylab="Density",main="")

lines(density(stimuli$logOLDW[stimuli$foreign == 1],cut=0),col = alpha("purple",0.6),lwd=2)
lines(density(stimuli$logOLDW[stimuli$foreign == 0],cut=0),col = alpha(words_col,1),lwd=2)
lines(density(stimuli$logOLDPW[!is.na(stimuli$logOLDPW)],cut=0),col = alpha(non_word_col[2],1),lwd=2)
lines(density(stimuli$logOLDRS[!is.na(stimuli$logOLDRS)],cut=0),col = alpha(non_word_col[1],1),lwd=2)
```

## Average results
```{r}
# Stage 5 = extra stage only present for pw
# Stage 6 = shared fifth stage in the manuscript
# Stage 7 = shared sixth stage in the manuscript

# Average stage duration per subject
subj_avg <- ddply(hsmm_dat,c("subjects_var","conditions"),summarize,
                  d1 = mean(durations_1),
                  d2 = mean(durations_2),
                  d3 = mean(durations_3),
                  d4 = mean(durations_4),
                  d5 = mean(durations_5),
                  d6 = mean(durations_6),
                  d7 = mean(durations_7))

# conditions averages + standard error over subjects
cond_avg <- ddply(subj_avg,c("conditions"),summarize,
                  sd1 = se(d1),
                  sd2 = se(d2),
                  sd3 = se(d3),
                  sd4 = se(d4),
                  sd5 = se(d5),
                  sd6 = se(d6),
                  sd7 = se(d7),
                  d1 = mean(d1),
                  d2 = mean(d2),
                  d3 = mean(d3),
                  d4 = mean(d4),
                  d5 = mean(d5),
                  d6 = mean(d6),
                  d7 = mean(d7))

cond_means <- t(as.matrix(cond_avg[,9:15]))
cond_ses <- t(as.matrix(cond_avg[,2:8]))
colnames(cond_means) <- cond_avg$conditions
colnames(cond_ses) <- cond_avg$conditions

# Average duration estimates (in ms)
cond_means

# Corresponding standard error
cond_ses

# Plots (topologies can be found in ./results/plots)
emptyPlot(c(1,7),
          c(50,400),
          xlab="Processing Stage",
          ylab="Stage Duration",
          bty="n",
          xaxt = "n",
          ymark=NULL,
          main="")

lines(1:7,cond_means[,1],col=non_word_col[2])
points(1:7,cond_means[,1],col=non_word_col[2],cex=1,pch=16)
errorBars(1:7,cond_means[,1],cond_ses[,1],lwd=1,col=non_word_col[2])

lines(c(1,2,3,4,6,7),cond_means[c(1,2,3,4,6,7),2],col=non_word_col[1])
points(c(1,2,3,4,6,7),cond_means[c(1,2,3,4,6,7),2],col=non_word_col[1],cex=1,pch=16)
errorBars(c(1,2,3,4,6,7),cond_means[c(1,2,3,4,6,7),2],cond_ses[c(1,2,3,4,6,7),2],lwd=1,col=non_word_col[1])

lines(c(1,2,3,4,6,7),cond_means[c(1,2,3,4,6,7),3],col=words_col)
points(c(1,2,3,4,6,7),cond_means[c(1,2,3,4,6,7),3],col=words_col,cex=1,pch=16)
errorBars(c(1,2,3,4,6,7),cond_means[c(1,2,3,4,6,7),3],cond_ses[c(1,2,3,4,6,7),3],lwd=1,col=words_col)

for(st in 1:7){
  text(st,max(cond_means[st,]),st,pos=3)
}

legend("topleft",
       c("Words","Random Strings","Pseudo-words"),
       col = c(words_col,non_word_col[1],non_word_col[2]),
       pch = 16,
       bty="n",
       xpd=NA,horiz=F)

bar_height <- 65
bar_spacing <- 25
emptyPlot(c(0,800),
          c(0,3*bar_height+3*bar_spacing),
          xlab="Time (in ms)",
          yaxt="n",
          bty="n",
          xmark = NULL,
          ymark=NULL,
          main="")

start_y <- 0
c_labels <- c("PW","RS","W")
for (c in c(2,1,3)){
  start_x <- 0
  for (s in 1:nrow(cond_means)){
    polygon(c(start_x:(start_x+cond_means[s,c]+1),
              rev(start_x:(start_x+cond_means[s,c]+1))),
            c(rep(start_y,length(start_x:(start_x+cond_means[s,c]+1))), 
              rev(rep(start_y+bar_height,length(start_x:(start_x+cond_means[s,c]+1))))),
            col=colrc[s],border = NA)
    
    start_x <- start_x + cond_means[s,c]
  }
  text(-10,start_y+(bar_height/2),c_labels[c],pos=2,xpd=NA)
  start_y <- start_y + bar_height + bar_spacing
}

# Check correlations between stage 5 and 6
cor(hsmm_dat$durations_65[hsmm_dat$conditions == "W"],
    hsmm_dat$durations_7[hsmm_dat$conditions == "W"])

cor(hsmm_dat$durations_65[hsmm_dat$conditions == "PW"],
    hsmm_dat$durations_7[hsmm_dat$conditions == "PW"])

cor(hsmm_dat$durations_65[hsmm_dat$conditions == "RS"],
    hsmm_dat$durations_7[hsmm_dat$conditions == "RS"])

```

## Trial-level results
```{r}
# Initial analyses revealed two clusters of words, so the conditions variable is
# adjusted to treat foreign words as a separate condition.
hsmm_dat$conditions <- as.character(hsmm_dat$conditions)
hsmm_dat$conditions[hsmm_dat$foreign == 1] <- "FW"
hsmm_dat$conditions <- as.factor(hsmm_dat$conditions)

# Per Stage models - estimated to maximize REML score
m1 <- bam(log(durations_1) ~ conditions + s(logGoogleFreq,by=conditions) +
                       s(trials_var,subjects_var,bs="fs",m=1) + 
                       s(logGoogleFreq,subjects_var,by=conditions,bs="fs",m=1),
          data=hsmm_dat)

m2 <- bam(log(durations_2) ~ conditions + s(logGoogleFreq,by=conditions) +
                       s(trials_var,subjects_var,bs="fs",m=1) + 
                       s(logGoogleFreq,subjects_var,by=conditions,bs="fs",m=1),
          data=hsmm_dat)

m3 <- bam(log(durations_3) ~ conditions + s(logGoogleFreq,by=conditions) +
                       s(trials_var,subjects_var,bs="fs",m=1) + 
                       s(logGoogleFreq,subjects_var,by=conditions,bs="fs",m=1),
          data=hsmm_dat)

m4 <- bam(log(durations_4) ~ conditions + s(logGoogleFreq,by=conditions) +
                       s(trials_var,subjects_var,bs="fs",m=1) + 
                       s(logGoogleFreq,subjects_var,by=conditions,bs="fs",m=1),
          data=hsmm_dat)

m5 <- bam(log(durations_65) ~ conditions + s(logGoogleFreq,by=conditions) +
                       s(trials_var,subjects_var,bs="fs",m=1) + 
                       s(logGoogleFreq,subjects_var,by=conditions,bs="fs",m=1),
          data=hsmm_dat)

m6 <- bam(log(durations_7) ~ conditions + s(logGoogleFreq,by=conditions) +
                       s(trials_var,subjects_var,bs="fs",m=1) + 
                       s(logGoogleFreq,subjects_var,by=conditions,bs="fs",m=1),
          data=hsmm_dat)

# Store fitted m6 values for later correlation analysis
m6_fitted <- m6$fitted.values

models <- list(m1,m2,m3,m4,m5,m6)

# Model criticism
for (mi in 1:length(models)) {
  m <- models[[mi]]

  # Residual plots do not suggest extreme violations of the
  # regression assumptions (residual_i ~ N(0,sigma); see Wood, 2017)
  par(mfrow=c(2,2))
  gam.check(m)
  par(mfrow=c(1,1))
  acf_resid(m)
}

# Analysis
par(mfrow=c(2,3))
for (mi in 1:length(models)) {
  m <- models[[mi]]
  
  report_stats(m)
  
  # Plots of the estimated effects of frequency evaluated over
  # the frequency range applicable to every word type.
  emptyPlot(c(0,25),
            c(min(fitted(m)),
              max(fitted(m))),
            xlab="Log(Google frequency)",
            ylab=paste0("Pred. Log(Stage ",mi," duration)"),
            bty="n")
  
  plot_smooth(m,
              view="logGoogleFreq",
              cond = list("conditions"="RS"),
              col=non_word_col[1],
              n.grid=500,
              xlim = c(min(stimuli$logFreqRS,na.rm = T),
                       max(stimuli$logFreqRS,na.rm = T)),
              add=T)
  
  plot_smooth(m,
              view="logGoogleFreq",
              cond = list("conditions"="PW"),
              col=non_word_col[2],
              n.grid=500,
              xlim = c(min(stimuli$logFreqPW,na.rm = T),
                       max(stimuli$logFreqPW,na.rm = T)),
              add=T)
  
  plot_smooth(m,
              view="logGoogleFreq",
              cond = list("conditions"="W"),
              col=words_col,
              n.grid=500,
              xlim = c(min(stimuli$logFreqW[stimuli$foreign == 0],na.rm = T),
                       max(stimuli$logFreqW[stimuli$foreign == 0],na.rm = T)),
              add=T)
  
  plot_smooth(m,
              view="logGoogleFreq",
              cond = list("conditions"="FW"),
              n.grid=500,
              col="purple",
              xlim = c(min(stimuli$logFreqW[stimuli$foreign == 1],na.rm = T),
                       max(stimuli$logFreqW[stimuli$foreign == 1],na.rm = T)),
              add=T)
}

# Difference certainty intervals
diffs <- list(c("W","FW"),
                c("W","PW"),
                c("W","RS"),
                c("FW","PW"),
                c("FW","RS"),
                c("PW","RS"))
  
mins <- c(max(min(hsmm_dat$logGoogleFreq[hsmm_dat$conditions == "W"]),
             min(hsmm_dat$logGoogleFreq[hsmm_dat$conditions == "FW"])),
          max(min(hsmm_dat$logGoogleFreq[hsmm_dat$conditions == "W"]),
             min(hsmm_dat$logGoogleFreq[hsmm_dat$conditions == "PW"])),
          max(min(hsmm_dat$logGoogleFreq[hsmm_dat$conditions == "W"]),
             min(hsmm_dat$logGoogleFreq[hsmm_dat$conditions == "RS"])),
          max(min(hsmm_dat$logGoogleFreq[hsmm_dat$conditions == "FW"]),
             min(hsmm_dat$logGoogleFreq[hsmm_dat$conditions == "PW"])),
          max(min(hsmm_dat$logGoogleFreq[hsmm_dat$conditions == "FW"]),
             min(hsmm_dat$logGoogleFreq[hsmm_dat$conditions == "RS"])),
          max(min(hsmm_dat$logGoogleFreq[hsmm_dat$conditions == "PW"]),
             min(hsmm_dat$logGoogleFreq[hsmm_dat$conditions == "RS"])))
  
maxs <- c(min(max(hsmm_dat$logGoogleFreq[hsmm_dat$conditions == "W"]),
             max(hsmm_dat$logGoogleFreq[hsmm_dat$conditions == "FW"])),
          min(max(hsmm_dat$logGoogleFreq[hsmm_dat$conditions == "W"]),
             max(hsmm_dat$logGoogleFreq[hsmm_dat$conditions == "PW"])),
          min(max(hsmm_dat$logGoogleFreq[hsmm_dat$conditions == "W"]),
             max(hsmm_dat$logGoogleFreq[hsmm_dat$conditions == "RS"])),
          min(max(hsmm_dat$logGoogleFreq[hsmm_dat$conditions == "FW"]),
             max(hsmm_dat$logGoogleFreq[hsmm_dat$conditions == "PW"])),
          min(max(hsmm_dat$logGoogleFreq[hsmm_dat$conditions == "FW"]),
             max(hsmm_dat$logGoogleFreq[hsmm_dat$conditions == "RS"])),
          min(max(hsmm_dat$logGoogleFreq[hsmm_dat$conditions == "PW"]),
             max(hsmm_dat$logGoogleFreq[hsmm_dat$conditions == "RS"])))

# Black line = estimated difference of the effect of frequency for two
# conditions. Only evaluated over the frequency range that is valid for both
# conditions.
# Shaded area: Bonferroni certainty interval around the difference (see Wood, 2017)
# Red bars: Frequency ranges in which the certainty interval does not include
# zero.
par(mfrow=c(2,3))
for (mi in 1:length(models)) {
  m <- models[[mi]]
  
  for (diff_i in 1:length(diffs)) {
    emptyPlot(c(0,25),
              c(-1,1),
              xlab="Log(Google frequency)",
              ylab=paste0("Est. Diff. in Pred. Log(Stage ",mi," duration)"),
              main=paste0(diffs[[diff_i]][1]," - ",diffs[[diff_i]][2]),
              bty="n")
    
    plot_diff(m,
              view="logGoogleFreq",
              comp = list(conditions=diffs[[diff_i]]),
              n.grid=500,
              xlim = c(mins[diff_i],
                       maxs[diff_i]),
              add=T,
              se=qnorm((1 - (0.05/36/2)))) # Bonferroni-corrected confidence intervals
  }
}

# Alternatively, we can also consider the AIC difference to check whether the
# three-way interaction offers a better account to the data:

# Reduced models - including only main effects of frequency
m1r <- bam(log(durations_1) ~ conditions + s(logGoogleFreq) +
                       s(trials_var,subjects_var,bs="fs",m=1) + 
                       s(logGoogleFreq,subjects_var,by=conditions,bs="fs",m=1),
          data=hsmm_dat)

m2r <- bam(log(durations_2) ~ conditions + s(logGoogleFreq) +
                       s(trials_var,subjects_var,bs="fs",m=1) + 
                       s(logGoogleFreq,subjects_var,by=conditions,bs="fs",m=1),
          data=hsmm_dat)

m3r <- bam(log(durations_3) ~ conditions + s(logGoogleFreq) +
                       s(trials_var,subjects_var,bs="fs",m=1) + 
                       s(logGoogleFreq,subjects_var,by=conditions,bs="fs",m=1),
          data=hsmm_dat)

m4r <- bam(log(durations_4) ~ conditions + s(logGoogleFreq) +
                       s(trials_var,subjects_var,bs="fs",m=1) + 
                       s(logGoogleFreq,subjects_var,by=conditions,bs="fs",m=1),
          data=hsmm_dat)

m5r <- bam(log(durations_65) ~ conditions + s(logGoogleFreq) +
                       s(trials_var,subjects_var,bs="fs",m=1) + 
                       s(logGoogleFreq,subjects_var,by=conditions,bs="fs",m=1),
          data=hsmm_dat)

m6r <- bam(log(durations_7) ~ conditions + s(logGoogleFreq) +
                       s(trials_var,subjects_var,bs="fs",m=1) + 
                       s(logGoogleFreq,subjects_var,by=conditions,bs="fs",m=1),
          data=hsmm_dat)

models_r <- list(m1r,m2r,m3r,m4r,m5r,m6r)

for(mi in 1:6){
  
  cat(paste0("AIC diff: ",AIC(models_r[[mi]]) - AIC(models[[mi]]),"\n"))
}

# Reset conditions to W, PW, and RS again:
hsmm_dat$conditions <- as.character(hsmm_dat$conditions)
hsmm_dat$conditions[hsmm_dat$conditions == "FW"] <- "W"
hsmm_dat$conditions <- as.factor(hsmm_dat$conditions)
```

### Shared model for all stages
```{r}
# Adjusted the condition variable again to treat foreign words as a separate condition.
hsmm_dat$conditions <- as.character(hsmm_dat$conditions)
hsmm_dat$conditions[hsmm_dat$foreign == 1] <- "FW"
hsmm_dat$conditions <- as.factor(hsmm_dat$conditions)

hsmm_dat_long <- NULL

for (stage in 1:6){
  if (stage == 5){
    # Take combined duration of shared and extra stage for pw
    stage_label <- "durations_65"
  } else if (stage == 6){
    stage_label <- "durations_7"
  } else {
    stage_label <- paste0("durations_",stage)
  }
  stage_hsmm_dat <- hsmm_dat
  stage_hsmm_dat$durations <- stage_hsmm_dat[,stage_label]
  stage_hsmm_dat$stage <- stage
  hsmm_dat_long <- rbind(hsmm_dat_long,stage_hsmm_dat)
}
hsmm_dat_long$stage_cond <- paste0(hsmm_dat_long$stage,"_",hsmm_dat_long$conditions)
hsmm_dat_long$stage_cond <- as.factor(hsmm_dat_long$stage_cond)
hsmm_dat_long$stage <- as.factor(hsmm_dat_long$stage)
hsmm_dat_long$conditions_ord <- as.ordered(hsmm_dat_long$conditions)
hsmm_dat_long <- droplevels(hsmm_dat_long)

if (T){
  m_shared <- bam(log(durations) ~ stage_cond + s(logGoogleFreq,by=stage_cond) +
                       s(trials_var,subjects_var,by=stage,bs="fs",m=1) + 
                       s(logGoogleFreq,subjects_var,by=stage_cond,bs="fs",m=1),
            data=hsmm_dat_long)
  
  m_shared_simpler <- bam(log(durations) ~ stage_cond + s(logGoogleFreq,by=stage)
                          + s(logGoogleFreq,by=conditions_ord) +
                       s(trials_var,subjects_var,by=stage,bs="fs",m=1) + 
                       s(logGoogleFreq,subjects_var,by=stage_cond,bs="fs",m=1),
            data=hsmm_dat_long)
  
  saveRDS(m_shared,file="./results/models/m_shared.RDS")
  saveRDS(m_shared_simpler,file="./results/models/m_shared_simpler.RDS")
} else {
  m_shared <- readRDS("./results/models/m_shared.RDS")
  m_shared_simpler <- readRDS("./results/models/m_shared_simpler.RDS")
}

par(mfrow=c(2,3))
for (mi in 1:6) {
  
  # Plots of the estimated effects of frequency evaluated over
  # the frequency range applicable to every word type.
  emptyPlot(c(0,25),
            c(min(fitted(m_shared)),
              max(fitted(m_shared))),
            xlab="Log(Google frequency)",
            ylab=paste0("Pred. Log(Stage ",mi," duration)"),
            bty="n")
  
  plot_smooth(m_shared,
              view="logGoogleFreq",
              cond = list("stage_cond"=paste0(mi,"_RS"),"stage"=mi),
              col=non_word_col[1],
              n.grid=500,
              xlim = c(min(stimuli$logFreqRS,na.rm = T),
                       max(stimuli$logFreqRS,na.rm = T)),
              add=T)
  
  plot_smooth(m_shared,
              view="logGoogleFreq",
              cond = list("stage_cond"=paste0(mi,"_PW"),"stage"=mi),
              col=non_word_col[2],
              n.grid=500,
              xlim = c(min(stimuli$logFreqPW,na.rm = T),
                       max(stimuli$logFreqPW,na.rm = T)),
              add=T)
  
  plot_smooth(m_shared,
              view="logGoogleFreq",
              cond = list("stage_cond"=paste0(mi,"_W"),"stage"=mi),
              col=words_col,
              n.grid=500,
              xlim = c(min(stimuli$logFreqW[stimuli$foreign == 0],na.rm = T),
                       max(stimuli$logFreqW[stimuli$foreign == 0],na.rm = T)),
              add=T)
  
  plot_smooth(m_shared,
              view="logGoogleFreq",
              cond = list("stage_cond"=paste0(mi,"_FW"),"stage"=mi),
              n.grid=500,
              col="purple",
              xlim = c(min(stimuli$logFreqW[stimuli$foreign == 1],na.rm = T),
                       max(stimuli$logFreqW[stimuli$foreign == 1],na.rm = T)),
              add=T)
}

# Check AIC
AIC(m_shared_simpler) - AIC(m_shared)

# Problems with this model: sigma should clearly be a function of stage.
plot(resid(m_shared),col=rep(colrc,each=nrow(hsmm_dat)))
plot(fitted(m_shared),resid(m_shared),col=rep(colrc,each=nrow(hsmm_dat)))
qqnorm(resid(m_shared))
qqline(resid(m_shared))

# Reset conditions to W, PW, and RS again:
hsmm_dat$conditions <- as.character(hsmm_dat$conditions)
hsmm_dat$conditions[hsmm_dat$conditions == "FW"] <- "W"
hsmm_dat$conditions <- as.factor(hsmm_dat$conditions)
```

## SUBTLEX trial-level
```{r}
# For the SUBTLEX analysis we fit the models to all words together
# since (as discussed in the manuscript) the SUBTLEX frequency does not suffer
# from inflated scores for the foreign words like the Google frequency measure.
# Otherwise the models are identical to the ones fitted above (of course no
# by=conditions).

# Per Stage models - estimated to maximize REML score
m1 <- bam(log(durations_1) ~ s(logSubtlex) +
                       s(trials_var,subjects_var,bs="fs",m=1) + 
                       s(logSubtlex,subjects_var,bs="fs",m=1),
          data=droplevels(hsmm_dat[hsmm_dat$conditions == "W",]))

m2 <- bam(log(durations_2) ~ s(logSubtlex) +
                       s(trials_var,subjects_var,bs="fs",m=1) + 
                       s(logSubtlex,subjects_var,bs="fs",m=1),
          data=droplevels(hsmm_dat[hsmm_dat$conditions == "W",]))

m3 <- bam(log(durations_3) ~ s(logSubtlex) +
                       s(trials_var,subjects_var,bs="fs",m=1) + 
                       s(logSubtlex,subjects_var,bs="fs",m=1),
          data=droplevels(hsmm_dat[hsmm_dat$conditions == "W",]))

m4 <- bam(log(durations_4) ~ s(logSubtlex) +
                       s(trials_var,subjects_var,bs="fs",m=1) + 
                       s(logSubtlex,subjects_var,bs="fs",m=1),
          data=droplevels(hsmm_dat[hsmm_dat$conditions == "W",]))

m5 <- bam(log(durations_65) ~ s(logSubtlex) +
                       s(trials_var,subjects_var,bs="fs",m=1) + 
                       s(logSubtlex,subjects_var,bs="fs",m=1),
          data=droplevels(hsmm_dat[hsmm_dat$conditions == "W",]))

m6 <- bam(log(durations_7) ~ s(logSubtlex) +
                       s(trials_var,subjects_var,bs="fs",m=1) + 
                       s(logSubtlex,subjects_var,bs="fs",m=1),
          data=droplevels(hsmm_dat[hsmm_dat$conditions == "W",]))

models <- list(m1,m2,m3,m4,m5,m6)

# Plot the estimated effects of SUBTLEX frequency
par(mfrow=c(2,3))
for (mi in 1:length(models)) {
  m <- models[[mi]]
  plot_smooth(m,
              view="logSubtlex",
              xlab="Log(SUBTLEX frequency)",
              ylab=paste0("Pred. Log(Stage ",mi," duration)"),
              col=words_col,
              n.grid=500)
}
```

## ERPs
```{r}
# CH 13 = Pz
# Ch 16 = Oz

s_rate <- 512
s_len <- 1000/512

par(mfrow=c(2,3))
for (ch in c(13,16)) {
  
  EEG_data <- matlab_dat$erp.data
  EEG_data <- EEG_data[ch,,]
  
  conds <- c("W","PW","RS")
  cond_cols <- c(words_col,non_word_col[2],non_word_col[1])
  
  for (cond_i in 3:1) {
    cond_col <- cond_cols[cond_i]
    cond <- conds[cond_i]
    
    # Compute critical values in condition-specific frequency distribution
    # corresponding to 4 quartiles.
    c_quantiles <- quantile(hsmm_dat$logGoogleFreq[hsmm_dat$conditions == cond],
                          probs=seq(0,1,(1/4)))
    
    # Early effects are of interest.
    emptyPlot(c(-100,300),
              c(-3,3),
              bty="n",
              xlab="Time (in ms)",
              ylab="Activation (mV)")
    
    # Average EEG signal per subject, time-point, and quantile.
    # Then compute average and standard error per time-point and
    # quantile and plot.
    cond_dat <- t(EEG_data[,matlab_dat$conditions == cond_i])
    cond_freq <- hsmm_dat$logGoogleFreq[matlab_dat$conditions == cond_i]
    cond_sub <- hsmm_dat$subjects_var[matlab_dat$conditions == cond_i]
    
    for (qi in 2:length(c_quantiles)) {
      
      if(qi == length(c_quantiles)) {
        qt_dat <- cond_dat[cond_freq <= c_quantiles[qi] &
                           cond_freq >= c_quantiles[(qi-1)],]
        
        qt_sub <- cond_sub[cond_freq <= c_quantiles[qi] &
                           cond_freq >= c_quantiles[(qi-1)]]
      } else {
        qt_dat <- cond_dat[cond_freq < c_quantiles[qi] &
                           cond_freq >= c_quantiles[(qi-1)],]
        qt_sub <- cond_sub[cond_freq < c_quantiles[qi] &
                           cond_freq >= c_quantiles[(qi-1)]]
      }
      
      all_sub_dat <- NULL
      
      for (sub in levels(hsmm_dat$subjects_var)){
        sub_dat <- qt_dat[qt_sub == sub,]
        aggr_sub_dat <- colMeans(sub_dat,na.rm=T)
        time <- seq(1,length(aggr_sub_dat)*s_len,by=s_len) - 200
        
        sub_dat <- data.frame("EEG"=aggr_sub_dat,
                              "time"=time,
                              "sub"=rep(sub,length(time)))
        
        all_sub_dat <- rbind(all_sub_dat,sub_dat)
      }
      
      aggr_all <- ddply(all_sub_dat,c("time"),summarize,
                        mEEG = mean(EEG),
                        seEEG = se(EEG))
      
      plot_error(aggr_all$time,
                 aggr_all$mEEG,
                 aggr_all$seEEG,
                 col=cond_col,
                 lty=(qi-1),
                 shade=T)
    }
  }
}
```

### Frequency Analysis
```{r}
# Adjusted the condition variable again to treat foreign words as a separate condition.
hsmm_dat$conditions <- as.character(hsmm_dat$conditions)
hsmm_dat$conditions[hsmm_dat$foreign == 1] <- "FW"
hsmm_dat$conditions <- as.factor(hsmm_dat$conditions)

# Stage 6 = Stage 5b for pseudo-words
m6 <- bam(log(durations_6) ~ conditions + s(logGoogleFreq,by=conditions) +
                       s(trials_var,subjects_var,bs="fs",m=1) + 
                       s(logGoogleFreq,subjects_var,by=conditions,bs="fs",m=1),
          data=hsmm_dat)

# Stage 65 = Stage 5a + 5b for pseudo-words
m65 <- bam(log(durations_65) ~ conditions + s(logGoogleFreq,by=conditions) +
                       s(trials_var,subjects_var,bs="fs",m=1) + 
                       s(logGoogleFreq,subjects_var,by=conditions,bs="fs",m=1),
           data=hsmm_dat)

# Stage 65a = Stage 5a for pseudo-words
hsmm_dat$durations_65a <- hsmm_dat$durations_6
hsmm_dat$durations_65a[hsmm_dat$condition == "PW"] <- hsmm_dat$durations_5[hsmm_dat$condition == "PW"]

m65a <- bam(log(durations_65a)  ~ conditions + s(logGoogleFreq,by=conditions) +
                       s(trials_var,subjects_var,bs="fs",m=1) + 
                       s(logGoogleFreq,subjects_var,by=conditions,bs="fs",m=1),
             data=hsmm_dat)

mRT <- bam(log(rts_var) ~ conditions + s(logGoogleFreq,by=conditions) +
                       s(trials_var,subjects_var,bs="fs",m=1) + 
                       s(logGoogleFreq,subjects_var,by=conditions,bs="fs",m=1),
           data=hsmm_dat)

# Plot RT effects
emptyPlot(xlim=c(0,25),ylim=c(6,7.1),
       xlab="Log(Google frequency)",
       ylab="Log(Reaction time)")
  
plot_smooth(mRT,
            view="logGoogleFreq",
            cond = list("conditions"="RS"),
            col=non_word_col[1],
            n.grid=500,
            xlim = c(min(stimuli$logFreqRS,na.rm = T),
                     max(stimuli$logFreqRS,na.rm = T)),
            add=T)

plot_smooth(mRT,
            view="logGoogleFreq",
            cond = list("conditions"="PW"),
            col=non_word_col[2],
            n.grid=500,
            xlim = c(min(stimuli$logFreqPW,na.rm = T),
                     max(stimuli$logFreqPW,na.rm = T)),
            add=T)

plot_smooth(mRT,
            view="logGoogleFreq",
            cond = list("conditions"="W"),
            col=words_col,
            n.grid=500,
            xlim = c(min(stimuli$logFreqW[stimuli$foreign == 0],na.rm = T),
                     max(stimuli$logFreqW[stimuli$foreign == 0],na.rm = T)),
            add=T)

plot_smooth(mRT,
            view="logGoogleFreq",
            cond = list("conditions"="FW"),
            n.grid=500,
            col=alpha("purple",0.6),
            xlim = c(min(stimuli$logFreqW[stimuli$foreign == 1],na.rm = T),
                     max(stimuli$logFreqW[stimuli$foreign == 1],na.rm = T)),
            add=T)

# Plot Stage 5 effects (5a + 5b for pseudo-words)
emptyPlot(xlim=c(0,25),ylim=c(4,7),
       xlab="Log(Google frequency)",
       ylab="Log(Stage duration)")
  
plot_smooth(m65,
            view="logGoogleFreq",
            cond = list("conditions"="RS"),
            col=non_word_col[1],
            n.grid=500,
            xlim = c(min(stimuli$logFreqRS,na.rm = T),
                     max(stimuli$logFreqRS,na.rm = T)),
            add=T)

plot_smooth(m65,
            view="logGoogleFreq",
            cond = list("conditions"="PW"),
            col=non_word_col[2],
            n.grid=500,
            xlim = c(min(stimuli$logFreqPW,na.rm = T),
                     max(stimuli$logFreqPW,na.rm = T)),
            add=T)

plot_smooth(m65,
            view="logGoogleFreq",
            cond = list("conditions"="W"),
            col=words_col,
            n.grid=500,
            xlim = c(min(stimuli$logFreqW[stimuli$foreign == 0],na.rm = T),
                     max(stimuli$logFreqW[stimuli$foreign == 0],na.rm = T)),
            add=T)

plot_smooth(m65,
            view="logGoogleFreq",
            cond = list("conditions"="FW"),
            n.grid=500,
            col=alpha("purple",0.6),
            xlim = c(min(stimuli$logFreqW[stimuli$foreign == 1],na.rm = T),
                     max(stimuli$logFreqW[stimuli$foreign == 1],na.rm = T)),
            add=T)

# Plot Stage 5 effects for 5a+5b, 5a, and 5b separately for pseudo-words
emptyPlot(xlim=c(0,25),ylim=c(4,7),
       xlab="Log(Google frequency)",
       ylab="Log(Stage duration)")

plot_smooth(m65,
            view="logGoogleFreq",
            cond = list("conditions"="PW"),
            col=non_word_col[2],
            n.grid=500,
            xlim = c(min(stimuli$logFreqPW,na.rm = T),
                     max(stimuli$logFreqPW,na.rm = T)),
            add=T)

plot_smooth(m6, # 5b
            view="logGoogleFreq",
            cond = list("conditions"="PW"),
            col=non_word_col[2],
            n.grid=500,
            xlim = c(min(stimuli$logFreqPW,na.rm = T),
                     max(stimuli$logFreqPW,na.rm = T)),
            add=T,lty=3)

plot_smooth(m65a,
            view="logGoogleFreq",
            cond = list("conditions"="PW"),
            col=non_word_col[2],
            n.grid=500,
            xlim = c(min(stimuli$logFreqPW,na.rm = T),
                     max(stimuli$logFreqPW,na.rm = T)),
            add=T,lty=2)

# Correlation between fixed effects for RTs and stage five, stage six, and sub-stage 5b
cor(m65$fitted,mRT$fitted)
cor(m6_fitted,mRT$fitted)
cor(m6$fitted,mRT$fitted)

# Reset conditions to W, PW, and RS again:
hsmm_dat$conditions <- as.character(hsmm_dat$conditions)
hsmm_dat$conditions[hsmm_dat$conditions == "FW"] <- "W"
hsmm_dat$conditions <- as.factor(hsmm_dat$conditions)
```