In pab2163/TOT_ERP: totERP Data Package

Mediation Analysis

Using Bayesian Multilevel Mediation

R Package by Matti Vuorre: https://github.com/mvuorre/bmlm

Assumed causal structure: tot states --> recall probability, mediated by ERP data

Load packages and data

library(bmlm)
library(plyr)
library(ggplot2)
library(gridExtra)
library(tidyverse)
library(qgraph)
library(bayesplot)

load("../data/long_clip.rda")

Reformatting for models

# Get Mean Amplitudes 250-700ms
timerange_clip <- filter(long_clip, timestamp >= 250) #filter only timestamps after 250ms


#Get average amplitude from 250ms on
avAMP_clip <- timerange_clip %>%
  group_by(id, chlabel, recall, tot, epoch) %>%
  summarise(value = mean(uv))


#filter only the electrodes we want for analysis
avAMP_pickELEC_clip <- filter(avAMP_clip, chlabel == "C1" | chlabel == "C2" | chlabel == "Cz" | chlabel == "CP1" | chlabel == "CP2" | chlabel == "CPz" | chlabel == "P1" | chlabel == "P2" | chlabel == "Pz")

#average across channels to get 1 measurement per trial
avAMP_pickELEC_avChan_clip <- avAMP_pickELEC_clip %>%
  group_by(id, recall, tot, epoch) %>%
  summarise(value = mean(value))

Recode / get data formatted to run bmlm models

# Dummy Code TOTs 
avAMP_pickELEC_avChan_clip$tot<- revalue(avAMP_pickELEC_avChan_clip$tot, c("yes"="1", "no"="0"))
head(avAMP_pickELEC_avChan_clip)


# Subject Center Variables 
iso <- isolate(d = avAMP_pickELEC_avChan_clip, by = "id",
                    value = "value")


# Necessary for Model to Work for Some Reason 
#write.csv(iso, file = "../data/singleTrialSubjectCentered.csv")
#iso <- read.csv("../data/singleTrialSubjectCentered.csv")

Fit the BMLM Model

iso_fit <- mlm(d = iso,
               id = "id",
               x = "tot",
               m = "value_cw",
               y = "recall", 
               binary_y = TRUE,
               iter = 10000, cores =4)

#Check MCMC
pars <- c("a", "b", "cp", "corrab", "me")
mcmc_trace(as.data.frame(iso_fit), pars = pars)


#Inspect Model

a2 <- mlm_summary(iso_fit)
head(a2)

Plots

mediation_chart <-  mlm_path_plot(iso_fit, level = .95, text = T,
              xlab = "TOT (Y/N)",
              mlab = "ERP Amplitude",
              ylab = "P(Recall)", digits = 2)

mlm_pars_plot(iso_fit, type = "violin", color = "dodgerblue2")


caterpillar <- mlm_pars_plot(iso_fit, pars = c("u_me","me"), type = "coef", level = .95,
              p_size = 3, p_shape = 18) + labs(title = "Subject-Specific Mediated Effects", x = 'Subject') +
  ylab('Mediated Effect')+ theme_bw() + 
  theme(axis.text.x=element_blank(), axis.title.y =element_text(size = 14), axis.title.x =element_text(size = 14),
        plot.title = element_text(size = 16, face = 'bold')) + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

caterpillar

Dev Version for Plotting

d <- iso  # Your data
d <- select(d, id, epoch, tot, recall, value, value_cw) %>% 
  rename(trial = epoch, amplitude = value, amplitude_cw = value_cw) %>% 
  arrange(id, trial) %>% 
  mutate(id = as.integer(as.factor(id)))

head(d)

# Some data figures
g1 <- d %>% 
  ggplot(aes(tot, amplitude_cw)) +
  geom_smooth(method = "lm", se = F) +
  geom_smooth(method = "lm", se = F, aes(group=id), size=.2, col = "black")

g2 <- d %>% 
  ggplot(aes(amplitude_cw, recall)) +
  geom_smooth(method = "glm", se = T, method.args = list(family=binomial())) +
  geom_smooth(method = "glm", se = F, method.args = list(family=binomial()),
              aes(group = id), size = .2, col = "black")
g1
g2


spaghetti <- mlm_spaghetti_plot(iso_fit, 
                                iso, 
                                id = "id", 
                                x = "tot", 
                                m = "amplitude_cw", 
                                y = "recall", 
                                # mx = "data",  # This is the new argument
                                binary_y = TRUE)


grid.arrange(spaghetti[[1]], g1, nrow = 1)
grid.arrange(spaghetti[[2]], g2, nrow = 1)