# this script is the same as analyse_03 but uses the ROI average within a post-hoc window

library(lme4)
library(brms)
library(ggdist)
library(dplyr)
library(purrr)
library(readr)
library(tidyr)
library(ggplot2)
library(scales)
library(patchwork)
library(ggnewscale)

ggplot2::theme_set(ggplot2::theme_classic() + theme(strip.background = element_rect(fill = "white")))

cong_cols <- c("#E69F00", "#009E73")

cong_cols_light <- sapply(cong_cols, function(x) {
  x_rgb <- as.numeric(col2rgb(x))
  x_li_rgb <- round(rowMeans(cbind(x_rgb*1.25, c(255, 255, 255)*0.75)))
  x_li <- rgb(x_li_rgb[1], x_li_rgb[2], x_li_rgb[3], maxColorValue=255)
  x_li
}, USE.NAMES=FALSE)

# function to normalise between 0 and 1
norm01 <- function(x, ...) (x-min(x, ...))/(max(x, ...)-min(x, ...))

# get the stimuli's percentage of name agreement values
stim <- read_csv("boss.csv", col_types = cols(perc_name_agree_denom_fq_inputs = col_number())) %>%
  select(filename, perc_name_agree_denom_fq_inputs) %>%
  rename(perc_name_agree = perc_name_agree_denom_fq_inputs)

# import the max electrode data from the preprocessing, and set up the variables for the model
d <- list.files("max_elec_data", pattern=".+\\.csv$", full.names=TRUE) %>%
  map_dfr(function(f) {
    read_csv(
      f,
      col_types=cols(
        date = col_date(format = "%d/%m/%Y"),
        trial_save_time = col_time(format = "%H:%M:%S"),
        subj_id = col_character(),
        stim_grp = col_integer(),
        resp_grp = col_integer(),
        sex = col_character(),
        trl_nr = col_integer(),
        item_nr = col_integer(),
        condition = col_character(),
        eeg_trigg = col_integer(),
        image = col_character(),
        string = col_character(),
        corr_ans = col_character(),
        resp = col_character(),
        acc = col_integer(),
        fix1_jitt_flip = col_integer(),
        fix2_jitt_flip = col_integer(),
        event = col_integer(),
        is_block_start = col_logical(),
        is_block_end = col_logical(),
        .default = col_double()
      )
    ) %>%
      select(
        subj_id, stim_grp, resp_grp, item_nr, condition,
        image, string, acc, rt, uV, uV_roi
      )
  }) %>%
  left_join(stim, by=c("image" = "filename")) %>%
  mutate(
    cong_dev = as.numeric(scale(ifelse(condition=="A2", 0, 1), center=TRUE, scale=FALSE)),
    cong_dum_incong = as.numeric(condition=="A1"),  # 0 is at incongruent
    cong_dum_cong = as.numeric(condition=="A2"),  # 0 is at congruent
    prop_agree = perc_name_agree/100,
    pred_norm = norm01(prop_agree, na.rm=TRUE),
    pred_norm_max = pred_norm - 1
  )

# fit the model ----
m <- lmer(
  uV_roi ~ cong_dev * pred_norm +
    (cong_dev * pred_norm | subj_id) +
    (cong_dev | image) +
    (1 | string),
  REML=FALSE,
  control = lmerControl(optimizer="bobyqa"),
  data=d
)

# the effect of the interaction
m_no_interact <- update(m, ~. -cong_dev:pred_norm)
eff_interact <- anova(m, m_no_interact)

# decompose interaction by congruency ----

# effect of predictability in incongruent trials
m_incong <- lmer(
  uV_roi ~ cong_dum_incong * pred_norm +
    (cong_dum_incong * pred_norm | subj_id) +
    (cong_dum_incong | image) +
    (1 | string),
  REML=FALSE,
  control = lmerControl(optimizer="bobyqa"),
  data=d
)

m_incong_no_pred <- update(m_incong, ~. -pred_norm)
eff_pred_incong <- anova(m_incong, m_incong_no_pred)

# effect of predictability in congruent trials
m_cong <- lmer(
  uV_roi ~ cong_dum_cong * pred_norm +
    (cong_dum_cong * pred_norm | subj_id) +
    (cong_dum_cong | image) +
    (1 | string),
  REML=FALSE,
  control = lmerControl(optimizer="bobyqa"),
  data=d
)

m_cong_no_pred <- update(m_cong, ~. -pred_norm)
eff_pred_cong <- anova(m_cong, m_cong_no_pred)

# decompose interaction by predictability ----

# effect of congruency at minimum level of predictability is in main model

# effect of congruency at maximum predictability
m_maxpred <- lmer(
  uV_roi ~ cong_dev * pred_norm_max +
    (cong_dev * pred_norm_max | subj_id) +
    (cong_dev | image) +
    (1 | string),
  REML=FALSE,
  control = lmerControl(optimizer="bobyqa"),
  data=d
)

m_maxpred_no_cong <- update(m_maxpred, ~. -cong_dev)
eff_cong_maxpred <- anova(m_maxpred, m_maxpred_no_cong)

# plot relationship ----

d_cells_pred <- tibble(
  pred_norm = rep(range(d$pred_norm), 2),
  perc_name_agree = rep(range(d$perc_name_agree), 2),
  cong_dev = rep(unique(d$cong_dev), each=2),
  condition = ifelse(cong_dev==min(cong_dev), "A2", "A1")
) %>%
  mutate(pred_uV = predict(m, newdata=., re.form=~0))

# get prediction intervals (no random slopes for feasibility)
m_no_rand_s <- lmer(
  uV_roi ~ cong_dev * pred_norm +
    (1 | subj_id) +
    (1 | image) +
    (1 | string),
  REML=FALSE,
  control = lmerControl(optimizer="bobyqa"),
  data=d
)

d_preds_boot <- expand_grid(
  cong_dev = unique(d$cong_dev),
  perc_name_agree = seq(min(d$perc_name_agree), max(d$perc_name_agree), 0.01)
) %>%
  mutate(
    prop_agree = perc_name_agree/100,
    pred_norm = norm01(prop_agree, na.rm=TRUE),
    condition = ifelse(cong_dev==min(cong_dev), "A2", "A1")
  )

boots <- bootMer(m_no_rand_s, nsim=5000, FUN=function(m_i) {
  predict(m_i, newdata=d_preds_boot, re.form=~0)
}, seed=3101, .progress="txt")

ci_norm <- function(x, level=c(0.025, 0.975)) {
  qnorm(p=level, mean=mean(x), sd=sd(x))
}

d_preds_boot <- d_preds_boot %>%
  mutate(
    pred_int_lwr = apply(boots$t, 2, ci_norm, level=.025),
    pred_int_upr = apply(boots$t, 2, ci_norm, level=.975)
  )

lmm_plot_scatter <- d %>%
  ggplot(aes(perc_name_agree, uV_roi, colour = condition)) +
  geom_point(shape=1, show.legend=FALSE, alpha=0.25) +
  scale_colour_manual(name = NA, labels = NULL, values = cong_cols_light) +
  new_scale_colour() +
  geom_line(aes(y=pred_uV, colour = condition), data=d_cells_pred, linewidth=1.75) +
  scale_colour_manual(name = "Picture-Word Congruency", labels = c("Congruent", "Incongruent"), values = cong_cols) +
  scale_y_continuous(breaks=scales::extended_breaks(n=6)) +
  labs(x = "Predictability (%)", y = "N1 Amplitude (µV)", tag="a")

lmm_plot_lines <- ggplot() +
  geom_ribbon(aes(x=perc_name_agree, ymin=pred_int_lwr, ymax=pred_int_upr, colour=condition), data=d_preds_boot, fill=NA, linetype="dashed", show.legend = FALSE) +
  geom_line(aes(x=perc_name_agree, y=pred_uV, colour=condition), data=d_cells_pred, linewidth=1.75) +
  scale_colour_manual(name = "Picture-Word Congruency", labels = c("Congruent", "Incongruent"), values = cong_cols) +
  scale_y_continuous(breaks=scales::extended_breaks(n=6)) +
  labs(x = "Predictability (%)", y = "N1 Amplitude (µV)", tag="b")

lmm_plot <- (lmm_plot_scatter | lmm_plot_lines) +
  plot_layout(guides = "collect") &
  theme(legend.position = "bottom", plot.background = element_blank()) &
  scale_x_continuous(expand=c(0, 0))

ggsave(file.path("figs", "05_lmm_summary_plot_roi.png"), lmm_plot, device="png", type="cairo", width=6.5, height=3.75, dpi=600)

ggsave(file.path("figs", "05_lmm_summary_plot_roi.pdf"), lmm_plot, device="pdf", width=6.5, height=3.75)

# evidences ratio for effect -----------------------------------------------

p <- c(
  set_prior("normal(-5, 10)", class="Intercept"),
  set_prior("normal(0, 5)", class="b", coef="cong_dev"),
  set_prior("normal(0, 5)", class="b", coef="pred_norm"),
  set_prior("normal(0, 5)", class="b", coef="cong_dev:pred_norm"),
  set_prior("lkj_corr_cholesky(1)", class="L")
)

m_b <- brm(
  uV_roi ~ cong_dev * pred_norm +
    (cong_dev * pred_norm | subj_id) +
    (cong_dev | image) +
    (1 | string),
  data = d,
  prior = p,
  chains = 5,
  cores = 5,
  iter = 5000,
  control = list(
    adapt_delta = 0.8,
    max_treedepth = 10
  ),
  refresh = 25,
  seed = 420
)

m_b_h <- hypothesis(m_b, "cong_dev:pred_norm>0")

m_b_h_hdi <- as_draws_df(m_b, variable="b_cong_dev:pred_norm", regex=FALSE) %>%
  rename(int_est = `b_cong_dev:pred_norm`) %>%
  median_hdi(int_est, .width=0.89)

m_b_pl <- as_draws_df(m_b, variable="b_cong_dev:pred_norm", regex=FALSE) %>%
  rename(int_est = `b_cong_dev:pred_norm`) %>%
  ggplot(aes(int_est, fill=after_stat(x>0))) +
  stat_slab(slab_colour="black", slab_size=0.4) +
  stat_pointinterval(
    .width=0.89, point_interval="median_hdi",
    point_size=1, size=1, show.legend = FALSE
  ) +
  geom_vline(xintercept=0, linetype="dashed") +
  geom_text(
    aes(x = x, y = y),
    data = tibble(x=-3, y=1),
    label = deparse(bquote(BF["01"]==.(round(1/m_b_h$hypothesis$Evid.Ratio, 2)))),
    parse = TRUE,
    hjust = 0,
    vjust = 1,
    size = 4.5
  ) +
  scale_y_continuous(expand = expansion(0, 0.04), limits = c(0, NA)) +
  coord_cartesian() +
  labs(
    x = "Congruency-Predictability Interaction (µV)",
    y = "Posterior Density"
  ) +
  theme(legend.position = "none") +
  scale_fill_manual(
    values = c("grey90", "#8B0000"),
    labels = c(bquote(H["0"]), bquote(H["1"]))
  )

ggsave(file.path("figs", "05_post.png"), width=3.5, height=1.75, device="png", type="cairo")
ggsave(file.path("figs", "05_post.pdf"), width=3.5, height=1.75, device="pdf")


lmm_plot_post <- wrap_plots(
  lmm_plot,
  (
    m_b_pl +
      labs(
        fill = NULL,
        tag = "c",
        x = "Congruency-Predictability\nInteraction"
      ) +
      theme(
        plot.background = element_blank(),
        legend.position = "bottom"
      )
  )
) +
  plot_layout(widths = c(2, 1))

ggsave(file.path("figs", "05_lmm_and_post.pdf"), lmm_plot_post, width=6.5, height=3, device="pdf")