AAT RTs baq, personal choice

No outliers, Hubert & Vandervieren (2008) method.

cat("\014")     # clean terminal



rm(list = ls()) # clean workspace
try(dev.off(), silent = TRUE) # close all plots
library(tidyverse)
library(readxl)
library(afex)
library(emmeans)
library(GGally)
library(robustbase)

my_dodge <- .5
exclude_bad_eeg <- TRUE
theme_set(theme_minimal())

a_posteriori <- function(afex_aov, sig_level = .05) {
  factors  <- as.list(rownames(afex_aov$anova_table))
  for (j in 1:length(factors)) {
    if (grepl(":", factors[[j]])) {
      factors[[j]] <- unlist(strsplit(factors[[j]], ":"))
    }
  }
  p_values <- afex_aov$anova_table$`Pr(>F)`
  for (i in 1:length(p_values)) {
    if (p_values[i] <= sig_level) {
      cat(rep("_", 100), '\n', sep = "")
      print(emmeans(afex_aov, factors[[i]], contr = "pairwise"))
    }
  }
}

eeg_check <- read_excel(file.path('..', '..', 'bad_channels_bqd.xlsx'))
eeg_check$num_id <- readr::parse_number(eeg_check$file)
bad_eeg  <- eeg_check$num_id[eeg_check$to_kill == 1]
RTs_name <- file.path('RTs_xtracted_from_EEG_personal_choice.csv')
RTs_data <- read.table(RTs_name, header = TRUE, strip.white = TRUE, sep = ",")
RTs_data <- RTs_data %>% separate(bin_name, c("cue","valence"), sep = "_", remove = FALSE)
RTs_data <- RTs_data %>% separate(set_name, c(NA, NA, "figure", "side"), sep = "_", remove = FALSE)
RTs_data$group[grepl("C1", RTs_data$set_name)] <- "Control"
RTs_data$group[grepl("C2", RTs_data$set_name)] <- "Mindfulness"
RTs_data$num_id  <- readr::parse_number(RTs_data$set_name)
RTs_data$num_id  <- factor(RTs_data$num_id)
RTs_data$log10_rt <- log10(RTs_data$rt)
RTs_data[sapply(RTs_data, is.character)] <- lapply(RTs_data[sapply(RTs_data, is.character)], as.factor)
if (exclude_bad_eeg) {
  RTs_data <- RTs_data[!(RTs_data$num_id %in% bad_eeg), ]
  }
lo_log10 <- adjbox(RTs_data$log10_rt, plot = FALSE)$fence[1]

The default of 'doScale' is FALSE now for stability;
  set options(mc_doScale_quiet=TRUE) to suppress this (once per session) message

up_log10 <- adjbox(RTs_data$log10_rt, plot = FALSE)$fence[2]
lo <- adjbox(RTs_data$rt, plot = FALSE)$fence[1]
up <- adjbox(RTs_data$rt, plot = FALSE)$fence[2]
RTs_data$rt_no_out <- RTs_data$rt
RTs_data$rt_no_out[RTs_data$rt < lo | RTs_data$rt > up] <- NA
RTs_data$rt_no_out_sec <- RTs_data$rt_no_out / 1000
RTs_data$log10_rt_no_out <- RTs_data$log10_rt
RTs_data$log10_rt_no_out[RTs_data$log10_rt < lo_log10 | RTs_data$log10_rt > up_log10] <- NA
RTs_data$cue <- fct_rev(RTs_data$cue)
write.csv(RTs_data,  file.path('RTs_xtracted_from_EEG_data_clean_personal_choice.csv'),  row.names = FALSE)

1 Participants

options(width = 100)
aggregate(data = RTs_data, num_id ~ group, function(num_id) length(unique(num_id)))

2 General description

options(width = 100)
summary(RTs_data)

              set_name      figure         side                      bin_name          cue       
 S19_C1_rect_right:  400   circ:11339   left :11755   Approach_Attractive:5721   Avoid   :11373  
 S02_C2_circ_left :  399   rect:11453   right:11037   Approach_Neutral   :5698   Approach:11419  
 S31_C1_rect_left :  399                              Avoid_Attractive   :5673                   
 S33_C1_circ_right:  399                              Avoid_Neutral      :5700                   
 S14_C2_circ_left :  398                                                                         
 S52_C2_rect_left :  398                                                                         
 (Other)          :20399                                                                         
       valence            rt                 group           num_id         log10_rt    
 Attractive:11394   Min.   : 301.8   Control    :11493   19     :  400   Min.   :2.480  
 Neutral   :11398   1st Qu.: 466.8   Mindfulness:11299   2      :  399   1st Qu.:2.669  
                    Median : 558.6                       31     :  399   Median :2.747  
                    Mean   : 616.2                       33     :  399   Mean   :2.764  
                    3rd Qu.: 687.5                       14     :  398   3rd Qu.:2.837  
                    Max.   :3265.6                       52     :  398   Max.   :3.514  
                                                         (Other):20399                  
   rt_no_out      rt_no_out_sec    log10_rt_no_out
 Min.   : 331.1   Min.   :0.3311   Min.   :2.502  
 1st Qu.: 470.7   1st Qu.:0.4707   1st Qu.:2.671  
 Median : 558.6   Median :0.5586   Median :2.747  
 Mean   : 598.4   Mean   :0.5984   Mean   :2.760  
 3rd Qu.: 681.6   3rd Qu.:0.6816   3rd Qu.:2.835  
 Max.   :1333.0   Max.   :1.3330   Max.   :3.178  
 NA's   :867      NA's   :867      NA's   :552    

3 Reaction times

3.1 All data

ggplot(
  RTs_data, aes(x = log10_rt, fill = cue, color = cue)) +
  geom_histogram(alpha = .3, position = "identity")

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

ggplot(
  RTs_data, aes(x = rt, fill = cue, color = cue)) +
  geom_histogram(alpha = .3, position = "identity")

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

3.2 No outliers

Hubert, M., & Vandervieren, E. (2008). An adjusted boxplot for skewed distributions. Computational Statistics & Data Analysis, 52(12), 5186–5201. https://doi.org/10.1016/J.CSDA.2007.11.008

ggplot(
  RTs_data, aes(x = log10_rt_no_out, fill = cue, color = cue)) +
  geom_histogram(alpha = .3, position = "identity")

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Warning: Removed 552 rows containing non-finite outside the scale range (`stat_bin()`).

ggplot(
  RTs_data, aes(x = rt_no_out, fill = cue, color = cue)) +
  geom_histogram(alpha = .3, position = "identity")

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Warning: Removed 867 rows containing non-finite outside the scale range (`stat_bin()`).

3.3 Chosen valences log10(RT)

options(width = 100)
valence_rep_anova <- aov_ez("num_id", "log10_rt_no_out", na.omit(RTs_data), within = c("cue", "valence"), between = c("group"))

Warning: More than one observation per design cell, aggregating data using `fun_aggregate = mean`.
To turn off this warning, pass `fun_aggregate = mean` explicitly.
Contrasts set to contr.sum for the following variables: group

cell_means <- valence_rep_anova$data$long |>
  group_by(group, valence, cue) |>
  summarise(log10_rt_no_out = mean(log10_rt_no_out))

`summarise()` has grouped output by 'group', 'valence'. You can override using the `.groups`
argument.

valenceafex_plot <-
  afex_plot(
    valence_rep_anova,
    x     = "valence",
    trace = "cue",
    panel = "group",
    error = "between",
    error_arg = list(width = .2, lwd = .75, col = 'gray30', alpha = .7),
    dodge = my_dodge,
    data_arg  = list(
      position = 
        position_jitterdodge(
          jitter.width  = .2, 
          # jitter.height = .1, 
          dodge.width   = my_dodge  ## needs to be same as dodge
        )),
    mapping = c("color"), data_alpha = .3,
    # point_arg = list(size = 2)
  )

Warning: Panel(s) show within-subjects factors, but not within-subjects error bars.
For within-subjects error bars use: error = "within"

nice(valence_rep_anova)

Anova Table (Type 3 tests)

Response: log10_rt_no_out
             Effect    df  MSE         F   ges p.value
1             group 1, 58 0.02      0.14  .002    .711
2               cue 1, 58 0.00 16.35 ***  .006   <.001
3         group:cue 1, 58 0.00      0.33 <.001    .568
4           valence 1, 58 0.00      0.02 <.001    .886
5     group:valence 1, 58 0.00      2.34 <.001    .131
6       cue:valence 1, 58 0.00      2.10 <.001    .152
7 group:cue:valence 1, 58 0.00      0.93 <.001    .339
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1

suppressWarnings(print(valenceafex_plot +
                         geom_point(data = cell_means,
                                    aes(x = valence, y = log10_rt_no_out, group = cue, color = cue),
                                    position = position_dodge(my_dodge), size = 2.5) +
                         facet_grid(~group)
))

a_posteriori(valence_rep_anova)

____________________________________________________________________________________________________
$emmeans
 cue      emmean      SE df lower.CL upper.CL
 Avoid      2.76 0.00840 58     2.75     2.78
 Approach   2.75 0.00848 58     2.74     2.77

Results are averaged over the levels of: group, valence 
Confidence level used: 0.95 

$contrasts
 contrast         estimate      SE df t.ratio p.value
 Avoid - Approach   0.0105 0.00259 58   4.044  0.0002

Results are averaged over the levels of: group, valence 

afex_plot(
  valence_rep_anova,
  x     = "valence",
  trace = "cue",
  panel = "group",
  error = "between",
  data_geom = geom_violin,
  error_arg = list(width = .2, lwd = .75, col = 'gray10', alpha = .8),
  dodge = my_dodge,
  line_arg = list(size = 0),
  mapping = c("color")
  ) + 
  geom_point(data = valence_rep_anova$data$long,
             aes(x = valence, y = log10_rt_no_out, group = cue, color = cue),
             position = position_jitterdodge(jitter.width = 0.2, dodge.width = my_dodge),
             alpha = .4) +
  geom_point(data = cell_means,
             aes(x = valence, y = log10_rt_no_out, group = cue, color = cue),
             position = position_dodge(my_dodge), size = 2.5) +
  facet_grid(~group)

Warning: Panel(s) show within-subjects factors, but not within-subjects error bars.
For within-subjects error bars use: error = "within"

3.4 Chosen valences RT

options(width = 100)
valence_rep_anova <- aov_ez("num_id", "rt_no_out", na.omit(RTs_data), within = c("cue", "valence"), between = c("group"))

Warning: More than one observation per design cell, aggregating data using `fun_aggregate = mean`.
To turn off this warning, pass `fun_aggregate = mean` explicitly.
Contrasts set to contr.sum for the following variables: group

cell_means <- valence_rep_anova$data$long |>
  group_by(group, valence, cue) |>
  summarise(rt_no_out = mean(rt_no_out))

`summarise()` has grouped output by 'group', 'valence'. You can override using the `.groups`
argument.

valenceafex_plot <-
  afex_plot(
    valence_rep_anova,
    x     = "valence",
    trace = "cue",
    panel = "group",
    error = "between",
    error_arg = list(width = .2, lwd = .75, col = 'gray30', alpha = .7),
    dodge = my_dodge,
    data_arg  = list(
      position = 
        position_jitterdodge(
          jitter.width  = .2, 
          # jitter.height = .1, 
          dodge.width   = my_dodge  ## needs to be same as dodge
        )),
    mapping = c("color"), data_alpha = .3,
    # point_arg = list(size = 2)
  )

Warning: Panel(s) show within-subjects factors, but not within-subjects error bars.
For within-subjects error bars use: error = "within"

nice(valence_rep_anova)

Anova Table (Type 3 tests)

Response: rt_no_out
             Effect    df      MSE         F   ges p.value
1             group 1, 58 34125.66      0.06 <.001    .815
2               cue 1, 58   807.12 14.08 ***  .006   <.001
3         group:cue 1, 58   807.12      0.25 <.001    .617
4           valence 1, 58   213.51      0.00 <.001    .993
5     group:valence 1, 58   213.51    3.84 + <.001    .055
6       cue:valence 1, 58   236.33      1.40 <.001    .241
7 group:cue:valence 1, 58   236.33      1.19 <.001    .280
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1

suppressWarnings(print(valenceafex_plot +
                         geom_point(data = cell_means,
                                    aes(x = valence, y = rt_no_out, group = cue, color = cue),
                                    position = position_dodge(my_dodge), size = 2.5) +
                         facet_grid(~group)
))

a_posteriori(valence_rep_anova)

____________________________________________________________________________________________________
$emmeans
 cue      emmean   SE df lower.CL upper.CL
 Avoid       604 12.2 58      580      629
 Approach    590 12.0 58      567      614

Results are averaged over the levels of: group, valence 
Confidence level used: 0.95 

$contrasts
 contrast         estimate   SE df t.ratio p.value
 Avoid - Approach     13.8 3.67 58   3.752  0.0004

Results are averaged over the levels of: group, valence 

afex_plot(
  valence_rep_anova,
  x     = "valence",
  trace = "cue",
  panel = "group",
  error = "between",
  data_geom = geom_violin,
  error_arg = list(width = .2, lwd = .75, col = 'gray10', alpha = .8),
  dodge = my_dodge,
  line_arg = list(size = 0),
  mapping = c("color")
  ) + 
  geom_point(data = valence_rep_anova$data$long,
             aes(x = valence, y = rt_no_out, group = cue, color = cue),
             position = position_jitterdodge(jitter.width = 0.2, dodge.width = my_dodge),
             alpha = .4) +
  geom_point(data = cell_means,
             aes(x = valence, y = rt_no_out, group = cue, color = cue),
             position = position_dodge(my_dodge), size = 2.5) +
  facet_grid(~group)

Warning: Panel(s) show within-subjects factors, but not within-subjects error bars.
For within-subjects error bars use: error = "within"

4 Generalized Linear Mixed-Effects Model

RTs in seconds because otherwise, for some weird reason, glmer() doesn’t work

options(width = 100)
library(lme4)
library(performance)
library(sjPlot)
library(report)
glmer_aat_full  <- glmer(rt_no_out_sec ~ valence * cue * group + (valence * cue | num_id), family = inverse.gaussian(log), data = RTs_data)
glmer_aat_noint <- glmer(rt_no_out_sec ~ valence + cue + group + (valence + cue | num_id), family = inverse.gaussian(log), data = RTs_data)

Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
  unable to evaluate scaled gradient
Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
  Model failed to converge: degenerate  Hessian with 1 negative eigenvalues

anova(glmer_aat_noint, glmer_aat_full)

Data: RTs_data
Models:
glmer_aat_noint: rt_no_out_sec ~ valence + cue + group + (valence + cue | num_id)
glmer_aat_full: rt_no_out_sec ~ valence * cue * group + (valence * cue | num_id)
                npar    AIC    BIC logLik deviance  Chisq Df Pr(>Chisq)  
glmer_aat_noint   11 -25116 -25028  12569   -25138                       
glmer_aat_full    19 -25116 -24964  12577   -25154 16.539  8    0.03528 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

summary(glmer_aat_full)

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: inverse.gaussian  ( log )
Formula: rt_no_out_sec ~ valence * cue * group + (valence * cue | num_id)
   Data: RTs_data

     AIC      BIC   logLik deviance df.resid 
-25116.4 -24964.5  12577.2 -25154.4    21906 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.9088 -0.6661 -0.1873  0.4424  8.5108 

Random effects:
 Groups   Name                       Variance  Std.Dev. Corr             
 num_id   (Intercept)                0.0045037 0.06711                   
          valenceNeutral             0.0003477 0.01865   0.03            
          cueApproach                0.0013680 0.03699  -0.25  0.27      
          valenceNeutral:cueApproach 0.0006846 0.02616   0.04 -0.50 -0.45
 Residual                            0.1122187 0.33499                   
Number of obs: 21925, groups:  num_id, 60

Fixed effects:
                                             Estimate Std. Error t value Pr(>|z|)    
(Intercept)                                 -0.504673   0.031217 -16.167   <2e-16 ***
valenceNeutral                               0.003688   0.008088   0.456   0.6483    
cueApproach                                 -0.026857   0.012486  -2.151   0.0315 *  
groupMindfulness                            -0.008168   0.044091  -0.185   0.8530    
valenceNeutral:cueApproach                   0.001933   0.011169   0.173   0.8626    
valenceNeutral:groupMindfulness             -0.016943   0.011446  -1.480   0.1388    
cueApproach:groupMindfulness                -0.001438   0.017666  -0.081   0.9351    
valenceNeutral:cueApproach:groupMindfulness  0.013699   0.015832   0.865   0.3869    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) vlncNt cApprc grpMnd vlnN:A vlnN:M cApp:M
valenceNtrl  0.082                                          
cueApproach -0.164  0.340                                   
grpMndflnss -0.708 -0.059  0.116                            
vlncNtrl:cA  0.032 -0.628 -0.505 -0.022                     
vlncNtrl:gM -0.059 -0.702 -0.239  0.081  0.440              
cApprch:grM  0.115 -0.239 -0.705 -0.165  0.355  0.341       
vlncNtr:A:M -0.022  0.440  0.354  0.033 -0.700 -0.628 -0.506

r2_nakagawa(glmer_aat_full)

# R2 for Mixed Models

  Conditional R2: 0.015
     Marginal R2: 0.001

# library(bayesplot)
# pp_check(glmer_aat_full)
check_predictions(glmer_aat_full)

Loading required namespace: statmod

check_model(glmer_aat_full)

plot_model(glmer_aat_full, type = "int")[[4]]

cue_aat <- emmeans(glmer_aat_full, "cue")

NOTE: Results may be misleading due to involvement in interactions

print(cue_aat)

 cue      emmean     SE  df asymp.LCL asymp.UCL
 Avoid    -0.511 0.0225 Inf    -0.555    -0.467
 Approach -0.534 0.0225 Inf    -0.578    -0.490

Results are averaged over the levels of: valence, group 
Results are given on the log (not the response) scale. 
Confidence level used: 0.95 

pairs(cue_aat)

 contrast         estimate      SE  df z.ratio p.value
 Avoid - Approach   0.0232 0.00764 Inf   3.034  0.0024

Results are averaged over the levels of: valence, group 
Results are given on the log (not the response) scale. 

plot(cue_aat, comparisons = TRUE)

# tab_model(glmer_aat_full)

5 This is the way, the Bayesian way…

library(bayesplot)

This is bayesplot version 1.11.1
- Online documentation and vignettes at mc-stan.org/bayesplot
- bayesplot theme set to bayesplot::theme_default()
   * Does _not_ affect other ggplot2 plots
   * See ?bayesplot_theme_set for details on theme setting

brm_aat_lognrm <- readRDS('brm_aat_lognrm.RData')
print(brm_aat_lognrm)

 Family: lognormal 
  Links: mu = identity; sigma = identity 
Formula: rt_no_out ~ group * valence * cue + (valence * cue | num_id) 
   Data: RTs_data (Number of observations: 21925) 
  Draws: 4 chains, each with iter = 10000; warmup = 5000; thin = 1;
         total post-warmup draws = 20000

Multilevel Hyperparameters:
~num_id (Number of levels: 60) 
                                            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)                                   0.15      0.02     0.13     0.19 1.00     2905     6482
sd(valenceNeutral)                              0.01      0.00     0.00     0.01 1.00     9670    10756
sd(cueAvoid)                                    0.04      0.01     0.03     0.05 1.00    12388    14727
sd(valenceNeutral:cueAvoid)                     0.01      0.01     0.00     0.02 1.00    10546     9656
cor(Intercept,valenceNeutral)                   0.18      0.41    -0.69     0.85 1.00    30160    14652
cor(Intercept,cueAvoid)                        -0.16      0.15    -0.45     0.14 1.00    14348    15092
cor(valenceNeutral,cueAvoid)                   -0.01      0.41    -0.78     0.77 1.00      700     2067
cor(Intercept,valenceNeutral:cueAvoid)          0.22      0.42    -0.69     0.87 1.00    28272    14860
cor(valenceNeutral,valenceNeutral:cueAvoid)    -0.09      0.45    -0.85     0.78 1.00    19843    16231
cor(cueAvoid,valenceNeutral:cueAvoid)          -0.08      0.42    -0.82     0.76 1.00    16121    14926

Regression Coefficients:
                                         Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept                                    6.34      0.03     6.28     6.40 1.00     1462     2477
groupMindfulness                            -0.01      0.04    -0.09     0.07 1.00     1466     2205
valenceNeutral                               0.01      0.01    -0.00     0.02 1.00    14580    16037
cueAvoid                                     0.03      0.01     0.01     0.05 1.00     9822    13545
groupMindfulness:valenceNeutral             -0.00      0.01    -0.02     0.01 1.00    15206    15945
groupMindfulness:cueAvoid                   -0.00      0.01    -0.03     0.03 1.00    10089    13650
valenceNeutral:cueAvoid                     -0.00      0.01    -0.02     0.01 1.00    12637    15387
groupMindfulness:valenceNeutral:cueAvoid    -0.01      0.01    -0.04     0.01 1.00    13591    15345

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.24      0.00     0.24     0.24 1.00    43399    14857

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

pp_check(brm_aat_lognrm, ndraws = 100)

pp_check(brm_aat_lognrm, type = "hist", ndraws = 19)

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(brm_aat_lognrm, type = "stat", stat = "mean", ndraws = 5000)

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(brm_aat_lognrm, type = "stat", stat = "min", ndraws = 5000)

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(brm_aat_lognrm, type = "stat", stat = "max", ndraws = 5000)

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.