Optimism_Worry

R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

summary(cars)

##      speed           dist       
##  Min.   : 4.0   Min.   :  2.00  
##  1st Qu.:12.0   1st Qu.: 26.00  
##  Median :15.0   Median : 36.00  
##  Mean   :15.4   Mean   : 42.98  
##  3rd Qu.:19.0   3rd Qu.: 56.00  
##  Max.   :25.0   Max.   :120.00

Including Plots

You can also embed plots, for example:

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.

title: “PerfectArtDiff” author: “Jutta Stahl” date: “September - 2024” output: html_document: toc: yes code_folding: hide toc_float: yes toc_depth: 3 number_sections: yes style: null pdf_document: toc: yes toc_depth: ‘3’ warning: False —

Code preparation and data preparation

Library & path

library(dplyr)
library(tidyr)
library(lme4)
library(lmerTest)
library(ggplot2)
library(afex)
library(apa)
library(emmeans)
library(knitr)
library(kableExtra)
library(ggsci)
library(cowplot)
library(texreg)
library(ggeffects)
library(interactions)
library(lsmeans)
library(flextable)
library(sjPlot)
library(scales)
library(Hmisc)
library(corrplot)
library(performance)
library(see)
library(PerformanceAnalytics)
library(tseries)
library(forcats)
library(readxl)
library(stringr)
library(parameters)
library(see)         
library(sjPlot)
library(interactions)
library(effectsize)

setwd(dir = "/Users/faezehsadipour/Desktop")

# Suppress summaries info
options(dplyr.summarise.inform = FALSE)
sessionInfo()

## R version 4.4.1 (2024-06-14)
## Platform: aarch64-apple-darwin20
## Running under: macOS Sonoma 14.0
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: Asia/Tehran
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] effectsize_1.0.0           parameters_0.26.0         
##  [3] stringr_1.5.1              readxl_1.4.3              
##  [5] forcats_1.0.0              tseries_0.10-58           
##  [7] PerformanceAnalytics_2.0.4 xts_0.14.0                
##  [9] zoo_1.8-12                 see_0.11.0                
## [11] performance_0.14.0         corrplot_0.94             
## [13] Hmisc_5.1-3                scales_1.3.0              
## [15] sjPlot_2.8.17              flextable_0.9.6           
## [17] lsmeans_2.30-0             interactions_1.2.0        
## [19] ggeffects_1.7.1            texreg_1.39.4             
## [21] cowplot_1.1.3              ggsci_3.2.0               
## [23] kableExtra_1.4.0           knitr_1.46                
## [25] emmeans_1.10.4             apa_0.3.4                 
## [27] afex_1.4-1                 ggplot2_3.5.1             
## [29] lmerTest_3.1-3             lme4_1.1-35.5             
## [31] Matrix_1.7-0               tidyr_1.3.1               
## [33] dplyr_1.1.4               
## 
## loaded via a namespace (and not attached):
##   [1] rstudioapi_0.16.0       jsonlite_1.8.8          datawizard_1.1.0       
##   [4] magrittr_2.0.3          estimability_1.5.1      nloptr_2.1.1           
##   [7] rmarkdown_2.27          ragg_1.3.2              vctrs_0.6.5            
##  [10] minqa_1.2.8             askpass_1.2.0           base64enc_0.1-3        
##  [13] htmltools_0.5.8.1       curl_5.2.1              broom_1.0.5            
##  [16] cellranger_1.1.0        Formula_1.2-5           sjmisc_2.8.10          
##  [19] TTR_0.24.4              sass_0.4.9              parallelly_1.38.0      
##  [22] bslib_0.8.0             htmlwidgets_1.6.4       plyr_1.8.9             
##  [25] cachem_1.1.0            uuid_1.2-1              lifecycle_1.0.4        
##  [28] pkgconfig_2.0.3         sjlabelled_1.2.0        R6_2.5.1               
##  [31] fastmap_1.2.0           future_1.34.0           digest_0.6.36          
##  [34] numDeriv_2016.8-1.1     colorspace_2.1-0        furrr_0.3.1            
##  [37] textshaping_0.4.0       fansi_1.0.6             httr_1.4.7             
##  [40] abind_1.4-8             compiler_4.4.1          fontquiver_0.2.1       
##  [43] withr_3.0.2             pander_0.6.5            htmlTable_2.4.3        
##  [46] backports_1.4.1         carData_3.0-5           highr_0.10             
##  [49] broom.mixed_0.2.9.5     MASS_7.3-60.2           openssl_2.2.0          
##  [52] sjstats_0.19.0          tools_4.4.1             foreign_0.8-86         
##  [55] quantmod_0.4.26         zip_2.3.1               nnet_7.3-19            
##  [58] glue_1.7.0              quadprog_1.5-8          nlme_3.1-164           
##  [61] grid_4.4.1              checkmate_2.3.2         cluster_2.1.6          
##  [64] reshape2_1.4.4          generics_0.1.3          gtable_0.3.5           
##  [67] MBESS_4.9.3             data.table_1.15.4       xml2_1.3.6             
##  [70] car_3.1-3               utf8_1.2.4              pillar_1.9.0           
##  [73] splines_4.4.1           lattice_0.22-6          tidyselect_1.2.1       
##  [76] fontLiberation_0.1.0    fontBitstreamVera_0.1.1 gridExtra_2.3          
##  [79] svglite_2.1.3           xfun_0.52               stringi_1.8.3          
##  [82] yaml_2.3.8              boot_1.3-30             evaluate_1.0.3         
##  [85] codetools_0.2-20        officer_0.6.6           gdtools_0.4.0          
##  [88] tibble_3.2.1            cli_3.6.2               rpart_4.1.23           
##  [91] xtable_1.8-4            systemfonts_1.1.0       munsell_0.5.1          
##  [94] jquerylib_0.1.4         Rcpp_1.0.13             globals_0.16.3         
##  [97] coda_0.19-4.1           parallel_4.4.1          bayestestR_0.16.0      
## [100] listenv_0.9.1           viridisLite_0.4.2       mvtnorm_1.3-1          
## [103] insight_1.3.0           purrr_1.0.2             rlang_1.1.3            
## [106] jtools_2.3.0

getwd()

## [1] "/Users/faezehsadipour/Desktop"

# Functions to plot the 2 and 3 way interactions with continous predictors
source("threeway_interactions.R")

## Customized functions
stderror <- function(x) sd(x)/sqrt(sum(!is.na(x)))

## Function to print lme result summary
tab_summary <- function(..., pred.labels, dv.labels, title = NULL) {
  tab_model(..., pred.labels = pred.labels, dv.labels = dv.labels, title = title, p.val = "satterthwaite",
            show.se = T, show.ci = 0.95, show.df =T, show.stat = T, show.obs = T, show.ngroups = T, 
            string.se = "SE", string.est = "Estimate", string.stat = "t", col.order = c("est", "std.error", "beta", "df.error", "stat", "p"))
  }

#1. Define bad participants (OK)
bad_participants <- c(4, 5, 13, 14, 25, 29, 48, 49, 75, 90, 96, 98)

# 2. Load questionnaire data (OK)
data.items <- read.csv2("/Users/faezehsadipour/Desktop/SE8ART/results_survey.csv")

# Define path and get filenames
path <- "/Users/faezehsadipour/Desktop/SE8ART/BEHAVIOUR"
files <- list.files(path = path, pattern = "SE8_HT[0-9]+\\.beh$", full.names = TRUE)

# Extract the numeric part from each file name
file_numbers <- as.numeric(sub("SE8_HT(\\d+)\\.beh$", "\\1", basename(files)))

# Optional: Clean invalid file names (just in case)
valid_indices <- !is.na(file_numbers)
files <- files[valid_indices]
file_numbers <- file_numbers[valid_indices]


data.orig <- bind_rows(lapply(seq_along(files), function(i) {
  file <- files[i]
  df <- read.table(file, header = TRUE)  # or use read.csv if that's more appropriate
  df$ID <- file_numbers[i]               # Add participant ID from filename
  df <- df %>%
    mutate(across(-ID, ~ as.numeric(.))) # Ensure numeric type, excluding ID
  return(df)
}))

# 4. Extract participant numbers
file_numbers <- as.numeric(sub("SE8_HT(\\d+)\\.beh$", "\\1", basename(files)))


# 6. Read and combine all files
data.orig <- bind_rows(lapply(seq_along(files), function(i) {
  data.orig <- data.orig %>%
  mutate(across(-ID, ~ as.numeric(.)))

  dplyr::last_dplyr_warnings()
data.orig <- data.orig %>%
  mutate(across(where(~ is.character(.) && !cur_column() %in% c("stimulus", "position")), 
                ~ as.numeric(.)))
sum(is.na(data.orig$RT))
cols_to_convert <- setdiff(names(data.orig), "ID")
data.orig[cols_to_convert] <- lapply(data.orig[cols_to_convert], function(x) as.numeric(x))
  # Read one file
  file_data <- read.csv(files[i], header = FALSE, sep = "\t")
  
  # Clean spaces
  file_data <- file_data %>% mutate_all(~ str_squish(.))
  
  # Add participant ID
  file_data$ID <- file_numbers[i]
  
  return(file_data)
}))

new_headers <- c(
  "block", "trial", "position", "stimulus", "ITI", 
  "force_1", "force_2", "force_3", "force_4", "force_5", "force_6", "force_7", "force_8",
  "RT_1", "RT_2", "RT_3", "RT_4", "RT_5", "RT_6", "RT_7", "RT_8",
  "TTP_1", "TTP_2", "TTP_3", "TTP_4", "TTP_5", "TTP_6", "TTP_7", "TTP_8",
  "firstresponse", "accuracy", "dummy", "RT", "RF", "TTP",
  "Det_force_1", "Det_force_2", "Det_force_3", "Det_force_4", "Det_force_5", "Det_force_6", "Det_force_7", "Det_force_8",
  "Det_RT_1", "Det_RT_2", "Det_RT_3", "Det_RT_4", "Det_RT_5", "Det_RT_6", "Det_RT_7", "Det_RT_8",
  "Det_TTP_1", "Det_TTP_2", "Det_TTP_3", "Det_TTP_4", "Det_TTP_5", "Det_TTP_6", "Det_TTP_7", "Det_TTP_8",
  "det_accu", "Det_RT", "DET_RF", "DET_TTP",
  "ID"
)

if (ncol(data.orig) == length(new_headers)) {
  colnames(data.orig) <- new_headers
}
data.orig <- data.orig %>%
  mutate(across(-ID, ~ as.numeric(.)))

length(unique(data.orig$ID))  

## [1] 91

if (ncol(data.orig) == length(new_headers)) {
  colnames(data.orig) <- new_headers
} else {
  stop("❌ Number of columns does not match expected structure.")
}

# === 6. Keep only blocks 9 to 15
data.orig <- data.orig %>%
 filter(block > 8 & block <= 15)

# Convert both ID columns to the same type
data.orig <- data.orig %>%
  mutate(ID = as.integer(ID))

data.items <- data.items %>%
  mutate(ID = as.integer(ID))

data.orig$ID <- as.integer(data.orig$ID)
data.items$ID <- as.integer(data.items$ID)

#adding data items to data orig 

data.orig <- left_join(data.orig, data.items, by = "ID")


#data_erp
path_erpagg <- "/Users/faezehsadipour/Desktop/SE8ART/Peaks.csv"
data.erpagg <- read.csv(path_erpagg, header = FALSE, sep = ";") %>%
  mutate_all(~ str_squish(.))
  new_headers <- c("ID", "Ne_Pcorrect_Selfeval", "Ne_Perror_Selfeval", 
                 "Pe_Perror_Selfeval", "Pe_Pcorrect_Selfeval")

if (ncol(data.erpagg) == length(new_headers)) {
  colnames(data.erpagg) <- new_headers
}
data.orig$ID <- as.integer(data.orig$ID)
data.erpagg$ID <- as.integer(data.erpagg$ID)
# View the cleaned data

#adding data erpagg to data orig 
data.orig <- left_join(data.orig, data.erpagg, by = "ID")
data.orig <- data.orig %>%
  relocate(ID, .before = everything())

Data housekeeping

data.orig <- data.orig %>%
  mutate(
    accuStr = case_when(
      accuracy == 1 ~ "correct",
      accuracy == 2 ~ "error",
      accuracy == 3 ~ "error_timeout",
      accuracy == 4 ~ "timeout",
      TRUE ~ NA_character_
    )
  )
data.orig <- data.orig %>%
  arrange(ID, trial) %>%  # ensure trial order
  mutate(
    Type_Nplus = lead(accuStr)  # next trial's accuracy string
  )
data.orig <- data.orig %>%
  arrange(ID, block, trial) %>%
  group_by(ID) %>%
  mutate(
    trial_index = row_number(),
    Type_Nplus = lead(accuStr)  # assign next trial's accuracy
  ) %>%
  ungroup()

data.orig <- data.orig %>% 
  mutate(
    # Binary flags for current trial
    error_b = ifelse(accuStr == "error", 1, 0),
    correct_b = ifelse(accuStr == "correct", 1, 0),
    error_timeout_b = ifelse(accuStr == "error_timeout", 1, 0),
    timeout_b = ifelse(accuStr == "timeout", 1, 0),
    too_slow_b = ifelse(accuStr %in% c("error_timeout", "timeout"), 1, 0),
    
    # Binary flags for next trial (n+1)
    error_nplus = ifelse(Type_Nplus == "error", 1, 0),
    correct_nplus = ifelse(Type_Nplus == "correct", 1, 0),
    error_timeout_nplus = ifelse(Type_Nplus == "error_timeout", 1, 0),
    timeout_nplus = ifelse(Type_Nplus == "timeout", 1, 0),
    too_slow_nplus = ifelse(Type_Nplus %in% c("error_timeout", "timeout"), 1, 0)
  )
  
# Clean data set: include only correct and error trials (exclude too slow & missing)
data.main.practice <- data.orig %>%
  filter(accuStr %in% c("correct", "error")) %>%
  mutate(
    accuStr = factor(accuStr, levels = c("error", "correct")),
    
  ) %>%
  ungroup()

data.main <- data.orig %>%

  filter(!RT %in% c("999999")) %>%           # Remove trials with invalid RTs
  filter(accuStr %in% c("error", "correct")) %>%   # Keep only trials with valid accuracy labels
  mutate(accuStr = factor(accuStr, levels = c("error", "correct"))) %>%          # Make  a factor (categorical variable)
  ungroup()


# Prepare for post-error analyses: exclude trials with no next RT
data.main.postbeh <- data.main.practice %>%
  arrange(ID, trial) %>%  # ensure trial order per subject
  mutate(RT_trial_nplus = lead(RT)) %>%
  filter(!RT_trial_nplus %in% c(999999)) %>%  # 999999 is numeric, not string
  mutate(
    accuStr = factor(accuStr),
     
  ) %>%
  ungroup()


questionnaire_worry <- c("msf.wX1.","msf.wX2.", "msf.wX3.", "msf.wX4.", "msf.wX5.", "msf.wX6.", "msf.wX7.", "msf.wX8.", "msf.wX9.", "msf.wX10.", "msf.wX11.", "msf.wX12.")
questionnaire_opt <- c("sop1.oX1.", "sop2.pX1.")
data.pers.cen <- data.items %>% 
  mutate(
    worry = rowMeans(across(all_of(questionnaire_worry)), na.rm = TRUE),
    opt = rowMeans(across(all_of(questionnaire_opt)), na.rm = TRUE),
worry_cen = worry - mean(worry, na.rm = TRUE),
 opt_cen = opt - mean(opt, na.rm = TRUE)
  )
data.pers.cen$ID <- as.character(data.pers.cen$ID)
data.main.postbeh$ID <- as.character(data.main.postbeh$ID)
data.main.postbeh <- data.main.postbeh %>%
  left_join(data.pers.cen %>% select(ID, opt, opt_cen, worry, worry_cen), by = "ID")
data.main$ID <- as.character(data.main$ID)
data.pers.cen$ID <- as.character(data.pers.cen$ID)
data.main <- data.main %>%
  left_join(data.pers.cen %>% select(ID, opt, opt_cen, worry, worry_cen), by = "ID")

#mean
colnames(data.items)

##  [1] "ID"        "age"       "gender"    "msf.wX1."  "msf.wX2."  "msf.wX3." 
##  [7] "msf.wX4."  "msf.wX5."  "msf.wX6."  "msf.wX7."  "msf.wX8."  "msf.wX9." 
## [13] "msf.wX10." "msf.wX11." "msf.wX12." "sop1.oX1." "sop2.pX1."

questionnaire_worry <- c("msf.wX1.","msf.wX2.", "msf.wX3.", "msf.wX4.", "msf.wX5.", "msf.wX6.", "msf.wX7.", "msf.wX8.", "msf.wX9.", "msf.wX10.", "msf.wX11.", "msf.wX12.")
questionnaire_opt <- c("sop1.oX1.", "sop2.pX1.")
data.pers.mean <- data.items %>% 
  mutate(
worry = rowMeans(across(all_of(questionnaire_worry)), na.rm = TRUE),
opt = rowMeans(across(all_of(questionnaire_opt)), na.rm = TRUE))

#Personality: only (non centered) 
selected_items <- data.items %>%
select(ID, age, gender, )
selected_cen <- data.pers.mean %>%
select(worry, opt) 
data.pers.all <- cbind(selected_items, selected_cen)
#SDT-parameter
data.orig <- data.orig %>%
  mutate(
    EU = ifelse(accuracy == 2 & det_accu == 5, "Error Undetected", NA),
    CU = ifelse(accuracy == 1 & det_accu == 6, "Correct Undetected", NA),
    ED = ifelse(accuracy == 2 & det_accu == 6, "Error Detected", NA),
    CD = ifelse(accuracy == 1 & det_accu == 5, "Correct Detected", NA)
  )

data.orig <- data.orig %>%
  mutate(
    signal = factor(case_when(
      !is.na(ED) ~ "ED",
      !is.na(EU) ~ "EU",
      !is.na(CU) ~ "CU",
      !is.na(CD) ~ "CD",
      TRUE ~ NA_character_
    )),
    hit = ifelse(!is.na(ED), 1, 0),
    fa = ifelse(!is.na(CU), 1, 0),
    miss = ifelse(!is.na(EU), 1, 0),
    cr = ifelse(!is.na(CD), 1, 0)
  ) %>%
  filter(!is.na(signal))

# Number of subjects
n_subj <- length(data.erpagg$ID)  

# Repeat subject-level variables twice (for correct and error conditions)
ID <- rep(data.erpagg$ID, 2)
worry_cen <- rep(data.pers.cen$worry_cen, 2)
opt_cen <- rep(data.pers.cen$opt_cen, 2)
worry <- rep(data.pers.all$worry, 2)
opt <- rep(data.pers.mean$opt, 2)

# Create accuracy factor (correct=1, error=2), length = 2*n_subj
acc <- factor(c(rep(1, n_subj), rep(2, n_subj)), labels = c("correct", "error"))
contrasts(acc) <- contr.sum(2) * (-0.5)  # sum contrast coding centered at 0

# Combine ERP amplitudes for Ne and Pe, concatenating correct and error values
Ne <- c(data.erpagg$Ne_Pcorrect_Selfeval, data.erpagg$Ne_Perror_Selfeval)
Pe <- c(data.erpagg$Pe_Pcorrect_Selfeval, data.erpagg$Pe_Perror_Selfeval)

# Create final data frame
matrix.mean.erp <- data.frame(
  ID, acc, worry_cen, opt_cen, worry, opt, Ne, Pe
) %>%
  mutate(
    Ne = as.numeric(Ne),
    Pe = as.numeric(Pe)
  )

# Define standard error function (handle NA properly)
stderror <- function(x) {
  sd(x, na.rm = TRUE) / sqrt(sum(!is.na(x)))
}

# Check that lengths match (should all be TRUE)
all_lengths_equal <- length(ID) == length(acc) && length(worry_cen) == length(acc) && length(Ne) == length(acc)
print(paste("All lengths equal:", all_lengths_equal))

## [1] "All lengths equal: TRUE"

matrix.mean.Ne3 <- matrix.mean.erp %>%
  group_by(ID, acc) %>%
  summarise(
    m.Ne = mean(Ne, na.rm = TRUE),
    se.Ne = stderror(Ne),
    worry_cen = first(worry_cen),
    opt_cen = first(opt_cen)
  ) %>%
  ungroup()

# 5. Check structure
str(matrix.mean.erp)

## 'data.frame':    182 obs. of  8 variables:
##  $ ID       : int  1 2 3 6 7 8 9 10 11 12 ...
##  $ acc      : Factor w/ 2 levels "correct","error": 1 1 1 1 1 1 1 1 1 1 ...
##   ..- attr(*, "contrasts")= num [1:2, 1] -0.5 0.5
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:2] "correct" "error"
##   .. .. ..$ : NULL
##  $ worry_cen: num  1.697 -0.886 -0.22 0.364 -0.22 ...
##  $ opt_cen  : num  -0.0934 -0.5934 -0.0934 -0.0934 -0.0934 ...
##  $ worry    : num  4.08 1.5 2.17 2.75 2.17 ...
##  $ opt      : num  4 3.5 4 4 4 4 4 3.5 3.5 4 ...
##  $ Ne       : num  -12.8 -14.2 -7.9 -5.6 -18.4 ...
##  $ Pe       : num  18.6 6.89 20.84 25.95 18.45 ...

# clean
rm(ID, opt , worry, opt_cen, worry_cen, acc, Ne, Pe)

Johnson-Neyman interval

(prepare data frame for predicted values)

# Convert `accuStr` to a factor with proper levels
data.main$accuStr <- factor(data.main$accuStr, levels = unique(data.main$accuStr))

# Check if levels are now correctly assigned
levels(data.main$accuStr)

## [1] "correct" "error"

limits <- max(c(abs(min(data.pers.cen$worry_cen)), abs(min(data.pers.cen$opt_cen)), max(data.pers.cen$worry_cen), max(data.pers.cen$opt_cen)))
b <- seq(-0.5-limits, 0.5+limits, 0.1)

# Limits are already defined earlier
data.pred1 <- data.frame(
  ID = rep(unique(data.main$ID), each = 2 * length(b) * length(b)),
  accuStr = rep(levels(data.main$accuStr), each = length(b) * length(b)),
  worry_cen = rep(b, each = length(b)),
  opt_cen = rep(b, length(b))
)

data.pred2 <- data.frame(
  ID = rep(unique(data.main$ID), each = 2 * length(b) * length(b)),
  accuStr = rep(levels(data.main$accuStr), each = length(b) * length(b)),
  worry_cen = rep(b, each = length(b)),
  opt_cen = rep(b, length(b))
)

# Merge and clean
data.pred <- rbind(data.pred1, data.pred2) %>%
  mutate(
    accuNum = ifelse(accuStr == "correct", -0.5, 0.5),
    accuStr = factor(accuStr)
  )

# Clean up temporary variables
rm(data.pred1, data.pred2)

jnp1 <- list()

Contrast coding

# Step 1: Convert to factor
data.main$accuStr <- factor(data.main$accuStr)

contrasts(data.main$accuStr) <- contr.sum(2) *(-0.5)

contrasts(data.main.postbeh$accuStr) <- contr.sum(2) *(-0.5)

Predictors

# full model
predictors.all <- c("(Intercept)", 
  "Response Type", 
  "Optimism", 
  "Worry", 
  "Response Type x Optimism", 
  "Evaluation Type x Optimism",  
  "Response Type x Worry", 
  "Worry x Optimism", 
  "Response Type x Worry x Optimism")

# Model with accuracy only
predictors.acc <- c(
  "(Intercept)", 
  "Response Type", 
  "Optimism", 
  "Worry", 
  "Response Type x Optimism", 
  "Response Type x Worry", 
  "Worry x Optimism", 
  "Response Type x Worry x Optimism"
)


# model only with personality variables
predictors.opt <- c("(Intercept)", 
  "Optimism", 
  "Worry", 
  "Worry x Optimism"
  
)

Demographics

Age

cat(
  "Age:\n",
  "Mean =", round(mean(data.pers.all$age, na.rm = TRUE), 2), 
  "+/-", round(sd(data.pers.all$age, na.rm = TRUE), 2), "\n",
  "Min =", round(min(data.pers.all$age, na.rm = TRUE), 2), "\n",
  "Max =", round(max(data.pers.all$age, na.rm = TRUE), 2)
)

## Age:
##  Mean = 22.68 +/- 5.68 
##  Min = 18 
##  Max = 59

Gender

paste("female:", sum(data.pers.all$gender == 1), "male:", sum(data.pers.all$gender == 2),"diverse:",sum(data.pers.all$gender == 0))

## [1] "female: 77 male: 14 diverse: 0"

Handedness?? I guess I don’t need it.

# Coefficient Oldfield
paste("right:", sum(data.erpagg$Latquo > 40), 
      "left:", sum(data.erpagg$Latquo < -40), 
    "ambidexter:", sum(data.erpagg$Latquo >= -40 & data.erpagg$Latquo <= 40))

## [1] "right: 0 left: 0 ambidexter: 0"

Scales and reliability

data.quest <- data.items %>% 
  mutate(
    worry = rowMeans(as.matrix(select(., msf.wX3., msf.wX4., msf.wX5., msf.wX6., 
                                      msf.wX7., msf.wX8., msf.wX9., msf.wX10., 
                                      msf.wX11., msf.wX12.)), na.rm = TRUE),

    opt = rowMeans(as.matrix(select(., sop1.oX1., sop2.pX1.)), na.rm = TRUE)
  ) %>% 
  ungroup()

# Optimism 
data.quest <- data.items %>% 
mutate(worry = mean(c(msf.wX3., msf.wX4., msf.wX5., msf.wX6., msf.wX7., msf.wX8., msf.wX9., msf.wX10., msf.wX11., msf.wX12. ), na.rm = T))%>%
mutate(opt = mean(c(sop1.oX1., sop2.pX1. ), na.rm = T)) %>% 
ungroup()

# Cronbachs alpha (conflicting - psych and scales, ggplot2) 
psych::alpha(subset(data.quest, select = c(msf.wX3., msf.wX4., msf.wX5., msf.wX6., msf.wX7., msf.wX8., msf.wX9., msf.wX10., msf.wX11., msf.wX12. )))

## 
## Reliability analysis   
## Call: psych::alpha(x = subset(data.quest, select = c(msf.wX3., msf.wX4., 
##     msf.wX5., msf.wX6., msf.wX7., msf.wX8., msf.wX9., msf.wX10., 
##     msf.wX11., msf.wX12.)))
## 
##   raw_alpha std.alpha G6(smc) average_r S/N  ase mean   sd median_r
##       0.93      0.94    0.94       0.6  15 0.01  2.4 0.95     0.61
## 
##     95% confidence boundaries 
##          lower alpha upper
## Feldt     0.91  0.93  0.95
## Duhachek  0.92  0.93  0.95
## 
##  Reliability if an item is dropped:
##           raw_alpha std.alpha G6(smc) average_r S/N alpha se  var.r med.r
## msf.wX3.       0.93      0.93    0.94      0.61  14    0.011 0.0087  0.62
## msf.wX4.       0.93      0.93    0.94      0.61  14    0.011 0.0078  0.62
## msf.wX5.       0.93      0.93    0.94      0.61  14    0.011 0.0088  0.60
## msf.wX6.       0.93      0.93    0.94      0.60  13    0.011 0.0089  0.60
## msf.wX7.       0.93      0.93    0.93      0.59  13    0.011 0.0088  0.60
## msf.wX8.       0.93      0.93    0.94      0.61  14    0.011 0.0092  0.62
## msf.wX9.       0.92      0.93    0.93      0.58  13    0.012 0.0086  0.59
## msf.wX10.      0.92      0.93    0.93      0.58  12    0.012 0.0081  0.59
## msf.wX11.      0.93      0.93    0.94      0.61  14    0.011 0.0080  0.62
## msf.wX12.      0.93      0.93    0.94      0.60  13    0.011 0.0089  0.60
## 
##  Item statistics 
##            n raw.r std.r r.cor r.drop mean   sd
## msf.wX3.  91  0.73  0.75  0.71   0.67  1.9 1.02
## msf.wX4.  91  0.77  0.77  0.74   0.71  3.3 1.12
## msf.wX5.  91  0.74  0.76  0.72   0.68  1.9 1.02
## msf.wX6.  91  0.81  0.80  0.78   0.75  2.5 1.32
## msf.wX7.  91  0.84  0.86  0.84   0.81  1.8 0.97
## msf.wX8.  91  0.77  0.76  0.72   0.70  2.9 1.42
## msf.wX9.  91  0.88  0.87  0.86   0.84  2.3 1.28
## msf.wX10. 91  0.88  0.88  0.87   0.84  2.6 1.30
## msf.wX11. 91  0.78  0.76  0.74   0.71  2.9 1.35
## msf.wX12. 91  0.79  0.80  0.77   0.74  2.0 1.06
## 
## Non missing response frequency for each item
##              1    2    3    4    5 miss
## msf.wX3.  0.49 0.22 0.20 0.09 0.00    0
## msf.wX4.  0.04 0.21 0.29 0.30 0.16    0
## msf.wX5.  0.46 0.29 0.14 0.11 0.00    0
## msf.wX6.  0.29 0.27 0.19 0.15 0.10    0
## msf.wX7.  0.48 0.27 0.16 0.08 0.00    0
## msf.wX8.  0.22 0.20 0.16 0.25 0.16    0
## msf.wX9.  0.37 0.26 0.13 0.18 0.05    0
## msf.wX10. 0.27 0.23 0.16 0.27 0.05    0
## msf.wX11. 0.20 0.24 0.16 0.26 0.13    0
## msf.wX12. 0.41 0.30 0.20 0.08 0.02    0

psych::alpha(subset(data.quest, select = c(sop1.oX1., sop2.pX1. )))

## Warning in psych::alpha(subset(data.quest, select = c(sop1.oX1., sop2.pX1.))): Some items were negatively correlated with the first principal component and probably 
## should be reversed.  
## To do this, run the function again with the 'check.keys=TRUE' option

## Some items ( sop1.oX1. ) were negatively correlated with the first principal component and 
## probably should be reversed.  
## To do this, run the function again with the 'check.keys=TRUE' option

## Warning in sqrt(Vtc): NaNs produced

## 
## Reliability analysis   
## Call: psych::alpha(x = subset(data.quest, select = c(sop1.oX1., sop2.pX1.)))
## 
##   raw_alpha std.alpha G6(smc) average_r   S/N  ase mean   sd median_r
##       -2.8      -2.9   -0.59     -0.59 -0.74 0.77  4.1 0.57    -0.59
## 
##     95% confidence boundaries 
##          lower alpha upper
## Feldt    -4.69 -2.76 -1.48
## Duhachek -4.26 -2.76 -1.25
## 
##  Reliability if an item is dropped:
##           raw_alpha std.alpha G6(smc) average_r   S/N alpha se var.r med.r
## sop1.oX1.     -0.49     -0.59    0.35     -0.59 -0.37       NA     0 -0.59
## sop2.pX1.     -0.70     -0.59    0.35     -0.59 -0.37       NA     0 -0.59
## 
##  Item statistics 
##            n raw.r std.r r.cor r.drop mean  sd
## sop1.oX1. 91  0.29  0.45   NaN  -0.59  5.1 1.1
## sop2.pX1. 91  0.60  0.45   NaN  -0.59  3.1 1.3
## 
## Non missing response frequency for each item
##              1    2    3    4    5    6   7 miss
## sop1.oX1. 0.00 0.03 0.03 0.19 0.41 0.24 0.1    0
## sop2.pX1. 0.04 0.41 0.20 0.18 0.11 0.07 0.0    0

Personality scores (mean and range)

paste("worry", round(mean(data.pers.all$worry),2), "+/-", round(sd(data.pers.all$worry), 2), "min:",round(min(data.pers.all$worry),2), "max:" , round(max(data.pers.all$worry),2))

## [1] "worry 2.39 +/- 0.94 min: 1 max: 4.5"

paste("opt", round(mean(data.pers.all$opt),2), "+/-", round(sd(data.pers.all$opt), 2),"min:",round(min(data.pers.all$opt),2), "max:" , round(max(data.pers.all$opt),2))

## [1] "opt 4.09 +/- 0.57 min: 2 max: 6"

Results Behaviour Trial n

Error rates

Descriptives

# matrix: ID, error rates, opt, worry
data.orig$ID <- as.character(data.orig$ID)
data.pers.cen$ID <- as.character(data.pers.cen$ID)
data.orig <- data.orig %>%
  mutate(ID = as.integer(as.character(ID)))
data.pers.cen <- data.pers.cen %>%
  mutate(ID = as.integer(as.character(ID)))
data.orig <- data.orig %>%
  left_join(data.pers.cen %>% select(ID, opt, opt_cen, worry, worry_cen), by = "ID")
names(data.orig)

##   [1] "ID"                   "block"                "trial"               
##   [4] "position"             "stimulus"             "ITI"                 
##   [7] "force_1"              "force_2"              "force_3"             
##  [10] "force_4"              "force_5"              "force_6"             
##  [13] "force_7"              "force_8"              "RT_1"                
##  [16] "RT_2"                 "RT_3"                 "RT_4"                
##  [19] "RT_5"                 "RT_6"                 "RT_7"                
##  [22] "RT_8"                 "TTP_1"                "TTP_2"               
##  [25] "TTP_3"                "TTP_4"                "TTP_5"               
##  [28] "TTP_6"                "TTP_7"                "TTP_8"               
##  [31] "firstresponse"        "accuracy"             "dummy"               
##  [34] "RT"                   "RF"                   "TTP"                 
##  [37] "Det_force_1"          "Det_force_2"          "Det_force_3"         
##  [40] "Det_force_4"          "Det_force_5"          "Det_force_6"         
##  [43] "Det_force_7"          "Det_force_8"          "Det_RT_1"            
##  [46] "Det_RT_2"             "Det_RT_3"             "Det_RT_4"            
##  [49] "Det_RT_5"             "Det_RT_6"             "Det_RT_7"            
##  [52] "Det_RT_8"             "Det_TTP_1"            "Det_TTP_2"           
##  [55] "Det_TTP_3"            "Det_TTP_4"            "Det_TTP_5"           
##  [58] "Det_TTP_6"            "Det_TTP_7"            "Det_TTP_8"           
##  [61] "det_accu"             "Det_RT"               "DET_RF"              
##  [64] "DET_TTP"              "age"                  "gender"              
##  [67] "msf.wX1."             "msf.wX2."             "msf.wX3."            
##  [70] "msf.wX4."             "msf.wX5."             "msf.wX6."            
##  [73] "msf.wX7."             "msf.wX8."             "msf.wX9."            
##  [76] "msf.wX10."            "msf.wX11."            "msf.wX12."           
##  [79] "sop1.oX1."            "sop2.pX1."            "Ne_Pcorrect_Selfeval"
##  [82] "Ne_Perror_Selfeval"   "Pe_Perror_Selfeval"   "Pe_Pcorrect_Selfeval"
##  [85] "accuStr"              "Type_Nplus"           "trial_index"         
##  [88] "error_b"              "correct_b"            "error_timeout_b"     
##  [91] "timeout_b"            "too_slow_b"           "error_nplus"         
##  [94] "correct_nplus"        "error_timeout_nplus"  "timeout_nplus"       
##  [97] "too_slow_nplus"       "EU"                   "CU"                  
## [100] "ED"                   "CD"                   "signal"              
## [103] "hit"                  "fa"                   "miss"                
## [106] "cr"                   "opt"                  "opt_cen"             
## [109] "worry"                "worry_cen"

# First: timeout summary
timeout_summary <- data.orig %>%
  group_by(ID) %>%
  summarise(
    total_trials = n(),
    timeout_abs = sum(RT == 999999.00, na.rm = TRUE),
    timeout_perc = round(100 * timeout_abs / total_trials, 2)
  )

sum(data.orig$RT == 999999.00, na.rm = TRUE) #to check for myself 

## [1] 0

matrix.mean.errorrate <- data.orig %>%
  group_by(ID) %>%
  summarise(
    error_abs = sum(accuStr == "error", na.rm = TRUE),
    correct_abs = sum(accuStr == "correct", na.rm = TRUE),
    total_trials = n(),
    error_perc = round(100 * error_abs / total_trials, 2),
    correct_perc = round(100 * correct_abs / total_trials, 2),
     too_slow_abs = sum(too_slow_b),
          too_slow_perc = mean(too_slow_b)*100,
          worry_cen = mean(worry_cen),
          opt_cen = mean(opt_cen)) %>%
ungroup()
 
#average.errorrate. <- matrix.mean.errorrate %>%
#group_by() %>%
#summarise(m.error_abs = mean(error_abs),
      
# average: error rates
##average.errorrate.EV <- matrix.mean.errorrate %>%
#group_by() %>%
#summarise(m.error_abs = mean(error_abs),
          #se.error_abs = stderror(error_abs),
          #m.error_perc = mean(error_perc),
          #se.error_perc = stderror(error_perc),
          #m.correct_abs = mean(correct_abs),
          #se.correct_abs = stderror(correct_abs),
          #m.correct_perc = mean(correct_perc),
          #se.correct_perc = stderror(correct_perc),
          #m.too_slow_abs = mean(too_slow_abs),
          #se.too_slow_abs = stderror(too_slow_abs),
          #m.too_slow_perc = mean(too_slow_perc),
          #se.too_slow_perc = stderror(too_slow_perc),
          #ps_cen = mean(ps_cen),
          #cm_cen = mean(cm_cen)
           # ) %>%
#ungroup()

# average: overall, error rates
se <- function(x) {
  sd(x, na.rm = TRUE) / sqrt(length(na.omit(x)))
}

average.errorrate.oa <- matrix.mean.errorrate %>%
group_by() %>%
summarise(m.error_abs = mean(error_abs),
          se.error_abs = stderror(error_abs),
          m.error_perc = mean(error_perc),
          se.error_perc = stderror(error_perc),
          m.correct_abs = mean(correct_abs),
          se.correct_abs = stderror(correct_abs),
          m.correct_perc = mean(correct_perc),
          se.correct_perc = stderror(correct_perc),
          m.too_slow_abs = mean(too_slow_abs),
          se.too_slow_abs = stderror(too_slow_abs),
          m.too_slow_perc = mean(too_slow_perc),
          se.too_slow_perc = stderror(too_slow_perc)
            ) %>%
ungroup()


paste("Error N:",round(average.errorrate.oa$m.error_abs,2), "+/-", round(average.errorrate.oa$se.error_abs,2), "Error rate: ", round(average.errorrate.oa$m.error_perc,2), "+/-", round(average.errorrate.oa$se.error_perc,2), "%")

## [1] "Error N: 27.37 +/- 1.73 Error rate:  7.75 +/- 0.54 %"

colnames(average.errorrate.oa)

##  [1] "m.error_abs"      "se.error_abs"     "m.error_perc"     "se.error_perc"   
##  [5] "m.correct_abs"    "se.correct_abs"   "m.correct_perc"   "se.correct_perc" 
##  [9] "m.too_slow_abs"   "se.too_slow_abs"  "m.too_slow_perc"  "se.too_slow_perc"

summary(average.errorrate.oa)

##   m.error_abs     se.error_abs    m.error_perc   se.error_perc   
##  Min.   :27.37   Min.   :1.733   Min.   :7.751   Min.   :0.5356  
##  1st Qu.:27.37   1st Qu.:1.733   1st Qu.:7.751   1st Qu.:0.5356  
##  Median :27.37   Median :1.733   Median :7.751   Median :0.5356  
##  Mean   :27.37   Mean   :1.733   Mean   :7.751   Mean   :0.5356  
##  3rd Qu.:27.37   3rd Qu.:1.733   3rd Qu.:7.751   3rd Qu.:0.5356  
##  Max.   :27.37   Max.   :1.733   Max.   :7.751   Max.   :0.5356  
##  m.correct_abs   se.correct_abs  m.correct_perc  se.correct_perc 
##  Min.   :338.2   Min.   :5.297   Min.   :92.25   Min.   :0.5356  
##  1st Qu.:338.2   1st Qu.:5.297   1st Qu.:92.25   1st Qu.:0.5356  
##  Median :338.2   Median :5.297   Median :92.25   Median :0.5356  
##  Mean   :338.2   Mean   :5.297   Mean   :92.25   Mean   :0.5356  
##  3rd Qu.:338.2   3rd Qu.:5.297   3rd Qu.:92.25   3rd Qu.:0.5356  
##  Max.   :338.2   Max.   :5.297   Max.   :92.25   Max.   :0.5356  
##  m.too_slow_abs se.too_slow_abs m.too_slow_perc se.too_slow_perc
##  Min.   :0      Min.   :0       Min.   :0       Min.   :0       
##  1st Qu.:0      1st Qu.:0       1st Qu.:0       1st Qu.:0       
##  Median :0      Median :0       Median :0       Median :0       
##  Mean   :0      Mean   :0       Mean   :0       Mean   :0       
##  3rd Qu.:0      3rd Qu.:0       3rd Qu.:0       3rd Qu.:0       
##  Max.   :0      Max.   :0       Max.   :0       Max.   :0

Error percentage

Error.lm <- lm(error_perc ~ opt_cen * worry_cen, data = matrix.mean.errorrate)
tab_model(Error.lm, dv.labels = "Error Rates [%]") 

	Error Rates [%]
Predictors	Estimates	CI	p
(Intercept)	7.67	6.53 – 8.81	<0.001
opt cen	0.19	-1.85 – 2.24	0.852
worry cen	0.25	-0.98 – 1.48	0.684
opt cen × worry cen	0.45	-1.56 – 2.47	0.655
Observations	91
R2 / R2 adjusted	0.006 / -0.029

Collinearity Error Percentage

check_collinearity(Error.lm)

## # Check for Multicollinearity
## 
## Low Correlation
## 
##               Term  VIF       VIF 95% CI adj. VIF Tolerance Tolerance 95% CI
##            opt_cen 1.14 [1.03,     1.80]     1.07      0.88     [0.56, 0.98]
##          worry_cen 1.14 [1.02,     1.81]     1.07      0.88     [0.55, 0.98]
##  opt_cen:worry_cen 1.00 [1.00,      Inf]     1.00      1.00     [0.00, 1.00]

# plot results
dum <- check_collinearity(Error.lm)
library(ggplot2)

ggplot(dum, aes(x = Term, y = VIF)) +
  geom_col() +
  geom_hline(yintercept = 5, linetype = "dashed", color = "red") +  # Warning line at VIF = 5
  theme_minimal() +
  labs(title = "Variance Inflation Factors (VIF)",
       y = "VIF Value",
       x = "Predictors")

Figures (Appendix)

# Fig.mean distributions   
  matrix.mean.errorrate$group <- "All"

ggplot(matrix.mean.errorrate, aes(x=group, y=error_perc)) +
  geom_violin(fill="#3df50a", color="#f20a19") +
  geom_boxplot(width=0.2) +
  geom_jitter(width=0.1) +
  labs(x = "Participants", y = "Error Rate [%]") +
  theme_classic()

### Statistical analyses

Signal detection theory parameter

Descriptives

# Base matrix with both predictors
matrix.mean.SDT <- data.orig %>%
  group_by(ID) %>%
  summarise(
    hit        = mean(hit),
    fa         = mean(fa),
    miss       = mean(miss),
    cr         = mean(cr),
    opt_cen    = mean(opt_cen),
    worry_cen  = mean(worry_cen)
  ) %>%
  ungroup() %>%
  mutate(
    dprime = qnorm(hit) - qnorm(fa),
    bias   = (-0.5) * (qnorm(hit) + qnorm(fa)),
    dprime = replace(dprime, is.infinite(dprime), NA),
    bias   = replace(bias, is.infinite(bias), NA)
  )

# 1. Main (Optimism-only)
matrix.mean.SDT.opt <- matrix.mean.SDT %>%
  select(ID, hit, fa, miss, cr, opt_cen, dprime, bias)

# 2. Exploratory (Optimism + Worry)
matrix.mean.SDT.exploratory <- matrix.mean.SDT


average.SDT.opt <- matrix.mean.SDT.opt %>%
  filter(!is.na(bias) & !is.na(dprime)) %>%
  summarise(
    m.bias    = mean(bias),
    se.bias   = stderror(bias),
    m.dprime  = mean(dprime),
    se.dprime = stderror(dprime)
  )

paste("sensitivity d′:", round(average.SDT.opt$m.dprime, 2), "+/-", round(average.SDT.opt$se.dprime, 2))

## [1] "sensitivity d′: 0.67 +/- 0.05"

paste("bias c:", round(average.SDT.opt$m.bias, 2), "+/-", round(average.SDT.opt$se.bias, 2))

## [1] "bias c: 1.94 +/- 0.04"

Certainty

Descriptives

# rating data

  matrix.mean.certain <-data.main.postbeh %>% 
  group_by(ID, accuStr, ) %>% 
  summarise(det_accu =  mean(det_accu),
            opt = mean(opt),
            worry = mean(worry)) %>%
  ungroup() 

# grand average : accuracy - RT
  average.certain <-  matrix.mean.certain %>%
  group_by(accuStr) %>%
  reframe(m.rating = mean(det_accu),
            se.rating = stderror(det_accu)) %>%
  ungroup()
  
paste("Rating: ",  average.certain$accuStr, round(average.certain$m.rating,2), "+/-", round(average.certain$se.rating,2))

## [1] "Rating:  error 5.34 +/- 0.05"   "Rating:  correct 4.52 +/- 0.05"

Statistical analyses

LME - SDT (bias’)

# linear model: bias, opt x worry 

bias.lm <- lm(bias ~ opt_cen * worry_cen, data = matrix.mean.SDT)


# 3. Define predictor labels (if needed)
predictors <- c("(Intercept)", "opt_cen", "worry_cen", "opt_cen:worry_cen")

ggplot(dum, aes(x = Term, y = VIF)) +
    geom_col() +
    geom_hline(yintercept = 5, linetype = "dashed", color = "red") +  # Warning line at VIF = 5
    theme_minimal() +
    labs(
      title = "Variance Inflation Factors (VIF)",
      y = "VIF Value",
      x = "Predictors"
    )

#### Collinearity bias

# no effect in LM

check_collinearity(bias.lm)

## # Check for Multicollinearity
## 
## Low Correlation
## 
##               Term  VIF       VIF 95% CI adj. VIF Tolerance Tolerance 95% CI
##            opt_cen 1.15 [1.03,     1.89]     1.07      0.87     [0.53, 0.97]
##          worry_cen 1.16 [1.03,     1.88]     1.08      0.87     [0.53, 0.97]
##  opt_cen:worry_cen 1.00 [1.00,      Inf]     1.00      1.00     [0.00, 1.00]

if (require("see")) {
  ggplot(dum, aes(x = Term, y = VIF)) +
    geom_col() +
    geom_hline(yintercept = 5, linetype = "dashed", color = "red") +  # Warning line at VIF = 5
    theme_minimal() +
    labs(
      title = "Variance Inflation Factors (VIF)",
      y = "VIF Value",
      x = "Predictors"
    )
}

LME - SDT (d’)

# linear model: d', opt x worry (only for Eval applicable) ?????
dprime.lm  <- lm(dprime ~ opt_cen * worry_cen, matrix.mean.SDT) 

tab_model(dprime.lm, dv.labels = "dprime [%]")

	dprime [%]
Predictors	Estimates	CI	p
(Intercept)	0.66	0.56 – 0.76	<0.001
opt cen	0.04	-0.13 – 0.21	0.667
worry cen	-0.01	-0.11 – 0.10	0.868
opt cen × worry cen	0.04	-0.13 – 0.21	0.625
Observations	76
R2 / R2 adjusted	0.006 / -0.036

Collinearity d prime

# no effect in LM

check_collinearity(dprime.lm)

## # Check for Multicollinearity
## 
## Low Correlation
## 
##               Term  VIF       VIF 95% CI adj. VIF Tolerance Tolerance 95% CI
##            opt_cen 1.15 [1.03,     1.89]     1.07      0.87     [0.53, 0.97]
##          worry_cen 1.16 [1.03,     1.88]     1.08      0.87     [0.53, 0.97]
##  opt_cen:worry_cen 1.00 [1.00,      Inf]     1.00      1.00     [0.00, 1.00]

# plot results
if (require("see")) {
  dum <- check_collinearity(dprime.lm )
  plot(dum)
}

Statistical analyses

LME - Certainty???

matrix.mean.certain <- data.main.postbeh %>%
  group_by(ID, accuStr) %>%
  summarise(
    det_accu = mean(det_accu, na.rm = TRUE),
    opt_cen = first(opt_cen),
    worry_cen = first(worry_cen)
  ) %>%
  ungroup()

Rating0.lme  <- lmer(det_accu ~ accuStr * opt_cen * worry_cen  + (1| ID), matrix.mean.certain)

tab_model(Rating0.lme, dv.labels = "Certainty [1..4] 1 = certain, 4 = uncertain")

	Certainty [1..4] 1 = certain, 4 = uncertain
Predictors	Estimates	CI	p
(Intercept)	4.93	4.85 – 5.01	<0.001
accuStr1	-0.82	-0.92 – -0.71	<0.001
opt cen	-0.11	-0.26 – 0.04	0.141
worry cen	0.09	0.00 – 0.18	0.049
accuStr1 × opt cen	0.14	-0.06 – 0.33	0.163
accuStr1 × worry cen	-0.07	-0.19 – 0.04	0.201
opt cen × worry cen	-0.01	-0.16 – 0.13	0.845
(accuStr1 × opt cen) × worry cen	0.02	-0.17 – 0.20	0.871
Random Effects
σ2	0.12
τ00 ID	0.08
ICC	0.41
N ID	91
Observations	182
Marginal R2 / Conditional R2	0.465 / 0.684

Collinearity Certainty????

check_collinearity(Rating0.lme)

## # Check for Multicollinearity
## 
## Low Correlation
## 
##                       Term  VIF       VIF 95% CI adj. VIF Tolerance
##                    accuStr 1.12 [1.03,     1.50]     1.06      0.89
##                    opt_cen 1.14 [1.04,     1.49]     1.07      0.88
##                  worry_cen 1.14 [1.04,     1.49]     1.07      0.88
##            accuStr:opt_cen 1.14 [1.04,     1.49]     1.07      0.88
##          accuStr:worry_cen 1.14 [1.04,     1.49]     1.07      0.88
##          opt_cen:worry_cen 1.00 [1.00, 9.71e+12]     1.00      1.00
##  accuStr:opt_cen:worry_cen 1.12 [1.03,     1.49]     1.06      0.89
##  Tolerance 95% CI
##      [0.67, 0.97]
##      [0.67, 0.96]
##      [0.67, 0.96]
##      [0.67, 0.96]
##      [0.67, 0.96]
##      [0.00, 1.00]
##      [0.67, 0.97]

# plot results
if (require("see")) {
  dum <- check_collinearity(Rating0.lme)
  plot(dum)
}

###Figures (Appendix)

# Fig.mean distributions   
  ggplot(matrix.mean.SDT, aes(x = 0, y=dprime)) +
  geom_violin(position=position_dodge(1)) +
  geom_point(position=position_dodge(1))+
  labs(x = "", y = "Sensitivity (d')") +
  geom_boxplot(position=position_dodge(1),width = 0.2)+
  theme_classic()

## Warning: Removed 15 rows containing non-finite outside the scale range
## (`stat_ydensity()`).

## Warning: Removed 15 rows containing non-finite outside the scale range
## (`stat_boxplot()`).

## Warning: Removed 15 rows containing missing values or values outside the scale range
## (`geom_point()`).

# Fig.mean distributions   
  ggplot(matrix.mean.SDT, aes(x = 0, y=bias)) +
  geom_violin(position=position_dodge(1)) +
  geom_point(position=position_dodge(1))+
  labs(x = "", y = "Bias (c)") +
  geom_boxplot(position=position_dodge(1),width = 0.2)+
  theme_classic()

## Warning: Removed 15 rows containing non-finite outside the scale range
## (`stat_ydensity()`).

## Warning: Removed 15 rows containing non-finite outside the scale range
## (`stat_boxplot()`).

## Warning: Removed 15 rows containing missing values or values outside the scale range
## (`geom_point()`).

Response Time

Descriptives

(Decision against median e.g. https://psyarxiv.com/3q5np/ )

# Summarise the data by ID and accuStr
matrix.mean.RT <- data.main %>%
  group_by(ID, accuStr ) %>%
  summarise(
    mRT = mean(RT, na.rm = TRUE),  # Calculate the mean of RT
    seRT = stderror(RT),           # Standard error of RT
    opt_cen = mean(opt_cen, na.rm = TRUE),
    worry_cen = mean(worry_cen, na.rm = TRUE)
  ) %>%
  ungroup()


average.RT.ID <- data.main %>%
  group_by(ID) %>%
  summarise(m.RT = mean(RT),
            se.RT = stderror(RT)) %>%
  ungroup()  
# grand average: experiment x accuracy - RT  
average.RT.AccEV <- matrix.mean.RT %>%
  group_by(accuStr, ) %>%
  summarise(m.RT = mean(mRT),
            se.RT.ae = sd(mRT)) %>%
  ungroup()
    
# grand average : accuracy - RT
average.RT.Acc <- matrix.mean.RT  %>%
  group_by(accuStr) %>%
  summarise(m.RT = mean(mRT),
            se.RT = stderror(mRT)) %>%
  ungroup()
  
# grand average : evaluation - RT  
#average.RT.EV <- matrix.mean.RT %>%
  #group_by() %>%
  #summarise(m.RT = mean(mRT),
           #se.RT = stderror(mRT)) %>%
  #ungroup()
  
average.RT.oa <- matrix.mean.RT  %>%
  group_by() %>%
  summarise(m.RT = mean(mRT),
            se.RT = stderror(mRT)) %>%
  ungroup()
    

paste("mean RT overall:", round(average.RT.oa$m.RT,2), "+/-",
round(average.RT.oa$se.RT,2), "ms")

## [1] "mean RT overall: 890.49 +/- 7.43 ms"

Statistical analyses

# Step 1: Set factor levels and contrasts
data.main.postbeh$accuStr <- factor(data.main.postbeh$accuStr, levels = c("error", "correct"))
contrasts(data.main.postbeh$accuStr) <- contr.sum(2) / 2  # sets: error = -0.5, correct = +0.5

# Step 2: Fit LME model with full 3-way interaction
RT_model1 <- lmer(RT ~ accuStr * opt_cen * worry_cen + (1 | ID), data = data.main.postbeh)

# Step 3: Simulated slope plot (3-way interaction)
RT_plot1 <- sim_slope_plot_3x(
  model = RT_model1,
  steps = 0.02,
  thresholding = TRUE
)

## Warning: Using an external vector in selections was deprecated in tidyselect 1.1.0.
## ℹ Please use `all_of()` or `any_of()` instead.
##   # Was:
##   data %>% select(contix)
## 
##   # Now:
##   data %>% select(all_of(contix))
## 
## See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: Using an external vector in selections was deprecated in tidyselect 1.1.0.
## ℹ Please use `all_of()` or `any_of()` instead.
##   # Was:
##   data %>% select(contiy)
## 
##   # Now:
##   data %>% select(all_of(contiy))
## 
## See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

# Step 4: Customize gradient fill and axes
RT_plot1 <- RT_plot1 +
  scale_fill_gradient2(
    limits = c(-150, 150),
    low = muted("blue"),
    mid = "white",
    high = muted("red"),
    name = "Accuracy Effect (error - correct)"
  ) +
  labs(
    title = "RT Interaction: Accuracy × Optimism × Worry",
    x = "Optimism (centered)",
    y = "Worry (centered)"
  ) +
  theme_classic() +
  theme(legend.position = "bottom")

## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.

# Step 5: Plot
cowplot::plot_grid(RT_plot1)

interact_plot(RT_model1, pred = "opt_cen", modx = "accuStr", mod2 = "worry_cen")

opt_range <- range(data.main.postbeh$opt_cen, na.rm = TRUE)
worry_range <- range(data.main.postbeh$worry_cen, na.rm = TRUE)
data.grid <- expand.grid(
  opt_cen = seq(opt_range[1], opt_range[2], length.out = 50),
  worry_cen = seq(worry_range[1], worry_range[2], length.out = 50),
  accuStr = factor(c("error", "correct"), levels = c("error", "correct"))
)

#### LME - RT #the table was not loaded.note for Faz  


# Accuracy x PSP x ECP
RT1.lme <- lmer(RT ~ accuStr * opt_cen * worry_cen + (accuStr | ID), data.main)

# binary - required for JN-procedure
 RT1.lmeb  <- lmer(RT  ~ accuStr * opt_cen * worry_cen + (accuStr | ID), data.main)
names(predictors.all) <- colnames(RT1.lme@vcov_beta)[1:length(predictors.all)]

# Table of all effects
names(predictors.all) <- colnames(RT1.lme@vcov_beta)

library(broom.mixed)
tidy(RT1.lme, effects = "fixed")

## # A tibble: 8 × 7
##   effect term                       estimate std.error statistic    df  p.value
##   <chr>  <chr>                         <dbl>     <dbl>     <dbl> <dbl>    <dbl>
## 1 fixed  (Intercept)                  884.        9.96   88.8     88.3 3.80e-88
## 2 fixed  accuStr1                      -2.33     10.2    -0.228   94.2 8.20e- 1
## 3 fixed  opt_cen                      -17.1      17.8    -0.957   88.1 3.41e- 1
## 4 fixed  worry_cen                     13.1      10.7     1.22    88.2 2.25e- 1
## 5 fixed  accuStr1:opt_cen              -1.33     18.3    -0.0726  92.3 9.42e- 1
## 6 fixed  accuStr1:worry_cen             5.33     11.0     0.486   93.6 6.28e- 1
## 7 fixed  opt_cen:worry_cen             33.6      17.7     1.90    89.8 6.05e- 2
## 8 fixed  accuStr1:opt_cen:worry_cen    19.7      18.4     1.07    98.0 2.87e- 1

tab_model(RT1.lme, 
          pred.labels = predictors.all, 
          dv.labels = "Response Time Effects")

	Response Time Effects
Predictors	Estimates	CI	p
(Intercept)	884.01	864.49 – 903.52	<0.001
Response Type	-2.33	-22.37 – 17.70	0.819
Optimism	-17.05	-51.99 – 17.89	0.339
Worry	13.06	-7.91 – 34.04	0.222
Response Type x Optimism	-1.33	-37.13 – 34.48	0.942
Evaluation Type x Optimism	5.33	-16.18 – 26.83	0.627
Response Type x Worry	33.62	-1.05 – 68.28	0.057
Worry x Optimism	19.67	-16.38 – 55.71	0.285
Random Effects
σ2	41380.23
τ00 ID	7588.01
τ11 ID.accuStr1	6582.06
ρ01 ID	0.39
ICC	0.14
N ID	91
Observations	37143
Marginal R2 / Conditional R2	0.007 / 0.149

tab_model(RT1.lme, 
          pred.labels = predictors.all, 
          dv.labels = "Response Time Effects")

	Response Time Effects
Predictors	Estimates	CI	p
(Intercept)	884.01	864.49 – 903.52	<0.001
Response Type	-2.33	-22.37 – 17.70	0.819
Optimism	-17.05	-51.99 – 17.89	0.339
Worry	13.06	-7.91 – 34.04	0.222
Response Type x Optimism	-1.33	-37.13 – 34.48	0.942
Evaluation Type x Optimism	5.33	-16.18 – 26.83	0.627
Response Type x Worry	33.62	-1.05 – 68.28	0.057
Worry x Optimism	19.67	-16.38 – 55.71	0.285
Random Effects
σ2	41380.23
τ00 ID	7588.01
τ11 ID.accuStr1	6582.06
ρ01 ID	0.39
ICC	0.14
N ID	91
Observations	37143
Marginal R2 / Conditional R2	0.007 / 0.149

  predictor.labels <- c(
  "(Intercept)" = "Baseline RT",
  "accuStr1" = "Accuracy: Error vs Correct",
  "opt_cen" = "Optimism (centered)",
  "worry_cen" = "Worry (centered)",
  "accuStr1:opt_cen" = "Accuracy × Optimism",
  "accuStr1:worry_cen" = "Accuracy × Worry",
  "opt_cen:worry_cen" = "Optimism × Worry",
  "accuStr1:opt_cen:worry_cen" = "3-Way Interaction"
)

Collinearity- RT

check_collinearity(RT1.lme )

## # Check for Multicollinearity
## 
## Low Correlation
## 
##                       Term  VIF   VIF 95% CI adj. VIF Tolerance
##                    accuStr 1.12 [1.11, 1.13]     1.06      0.89
##                    opt_cen 1.41 [1.39, 1.43]     1.19      0.71
##                  worry_cen 1.40 [1.38, 1.42]     1.18      0.71
##            accuStr:opt_cen 1.41 [1.39, 1.43]     1.19      0.71
##          accuStr:worry_cen 1.40 [1.38, 1.42]     1.18      0.71
##          opt_cen:worry_cen 1.25 [1.23, 1.26]     1.12      0.80
##  accuStr:opt_cen:worry_cen 1.38 [1.36, 1.39]     1.17      0.73
##  Tolerance 95% CI
##      [0.88, 0.90]
##      [0.70, 0.72]
##      [0.70, 0.72]
##      [0.70, 0.72]
##      [0.70, 0.72]
##      [0.79, 0.81]
##      [0.72, 0.74]

# plot results
if (require("see")) {
  x <- check_collinearity(RT1.lme)
  plot(x)
}

#### Post hoc - Acc

### Post hoc pairwise comparisons – Accuracy only
emmeans(RT1.lme, specs = "accuStr") %>% pairs()

## Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
## To enable adjustments, add the argument 'pbkrtest.limit = 37143' (or larger)
## [or, globally, 'set emm_options(pbkrtest.limit = 37143)' or larger];
## but be warned that this may result in large computation time and memory use.

## Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
## To enable adjustments, add the argument 'lmerTest.limit = 37143' (or larger)
## [or, globally, 'set emm_options(lmerTest.limit = 37143)' or larger];
## but be warned that this may result in large computation time and memory use.

## NOTE: Results may be misleading due to involvement in interactions

##  contrast        estimate   SE  df z.ratio p.value
##  correct - error     2.38 10.2 Inf   0.233  0.8158
## 
## Degrees-of-freedom method: asymptotic

emmeans(RT1.lme, specs = "accuStr", lmerTest.limit = 50000) %>% pairs()

## Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
## To enable adjustments, add the argument 'pbkrtest.limit = 37143' (or larger)
## [or, globally, 'set emm_options(pbkrtest.limit = 37143)' or larger];
## but be warned that this may result in large computation time and memory use.
## NOTE: Results may be misleading due to involvement in interactions

##  contrast        estimate   SE   df t.ratio p.value
##  correct - error     2.38 10.2 94.2   0.233  0.8163
## 
## Degrees-of-freedom method: satterthwaite

#### Johnson-Neyman: opt_cen × Accuracy

# Map opt_cen back to post-behavioral dataset
data.main.postbeh$opt_cen <- data.main$opt_cen[match(data.main.postbeh$ID, data.main$ID)]

# Contrast-code accuracy manually (bypasses parsing issues)
data.main.postbeh$accuNum <- ifelse(data.main.postbeh$accuStr == "error", 0.5, -0.5)

# Refit model using numeric-coded accuracy
RT1.simple.num <- lmer(RT ~ accuNum * opt_cen + (1 | ID), data = data.main.postbeh)

# Generate predicted RTs (optional)
data.pred <- data.main.postbeh %>%
  mutate(RTP = predict(RT1.simple.num, newdata = data.main.postbeh, allow.new.levels = TRUE))

# Johnson-Neyman interval (numeric output and plot)
RT1.jn <- interactions::johnson_neyman(
  model = RT1.simple.num,
  pred = "accuNum",      # contrast-coded: -0.5 = correct, 0.5 = error
  modx = "opt_cen",      # moderator = optimism
  plot = TRUE
)

# Summary
summary(RT1.simple.num)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: RT ~ accuNum * opt_cen + (1 | ID)
##    Data: data.main.postbeh
## 
## REML criterion at convergence: 501089.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.1488 -0.7296 -0.1222  0.6702  3.1894 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept)  6667     81.65  
##  Residual             41917    204.74  
## Number of obs: 37143, groups:  ID, 91
## 
## Fixed effects:
##                  Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)       888.541      8.796    96.016 101.017   <2e-16 ***
## accuNum            -1.379      4.082 37088.943  -0.338    0.735    
## opt_cen           -13.070     15.564    95.309  -0.840    0.403    
## accuNum:opt_cen    -3.069      6.971 37084.488  -0.440    0.660    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) accuNm opt_cn
## accuNum      0.196              
## opt_cen     -0.002 -0.008       
## accNm:pt_cn -0.008 -0.034  0.187

print(RT1.jn)

## JOHNSON-NEYMAN INTERVAL
## 
## The Johnson-Neyman interval could not be found. Is the p value for your
## interaction term below the specified alpha?

### Plot observed RT across optimism and accuracy (original factor-coded accuracy)
data.main.postbeh$accuStr <- factor(data.main.postbeh$accuStr, levels = c("correct", "error"))

plot_data <- data.main.postbeh %>%
  group_by(accuStr, opt_cen) %>%
  summarise(mean_RT = mean(RT, na.rm = TRUE), .groups = "drop")

ggplot(plot_data, aes(x = opt_cen, y = mean_RT, colour = accuStr, linetype = accuStr)) +
  geom_line(linewidth = 1.2) +
  theme_minimal() +
  labs(
    title = "Observed RT across Optimism and Accuracy",
    x = "Optimism (centered)",
    y = "Observed RT (ms)",
    colour = "Accuracy",
    linetype = "Accuracy"
  )

### Johnson-Neyman-style shaded plot
yrange <- range(data.main.postbeh$RT, na.rm = TRUE)

jnp1$RT.opt <- data.main.postbeh %>%
  group_by(accuStr, opt_cen, ID) %>%
  summarise(RT = mean(RT), .groups = "drop") %>%
  group_by(accuStr, opt_cen) %>%
  summarise(RT = mean(RT), .groups = "drop") %>%
  ggplot(aes(x = opt_cen, y = RT, colour = accuStr, linetype = accuStr)) +
  geom_rect(
    alpha = 0.2, fill = "lightgray", inherit.aes = FALSE,
    ymin = -Inf, ymax = Inf,
    xmin = min(data.main.postbeh$opt_cen),
    xmax = min(RT1.jn$bounds["Lower"])
  ) +
  geom_point(
    data = data.main.postbeh %>% select(ID, opt_cen) %>% distinct(),
    mapping = aes(x = opt_cen, y = yrange[1] + 0.01 * (yrange[2] - yrange[1])),
    inherit.aes = FALSE, alpha = 0.2,
    position = position_jitter(height = 0.01 * (yrange[2] - yrange[1]))
  ) +
  geom_line(linewidth = 1.1) +
  labs(
    x = "Optimism (centered)",
    y = "Observed Response Time [ms]",
    colour = "Accuracy",
    linetype = "Accuracy"
  ) +
  scale_color_grey(start = 0, end = 0.2) +
  scale_linetype_manual(values = c("solid", "longdash")) +
  theme_classic()

# Display final plot
jnp1$RT.opt

# Get observed RT per condition and optimism level
plot_data <- data.main.postbeh %>%
  group_by(accuStr, opt_cen) %>%
  summarise(mean_RT = mean(RT, na.rm = TRUE), .groups = "drop")

# Johnson-Neyman bounds from RT1.jn
jn_lower <- RT1.jn$bounds["Lower"]
jn_upper <- RT1.jn$bounds["Upper"]

# Base plot
ggplot(plot_data, aes(x = opt_cen, y = mean_RT, colour = accuStr, linetype = accuStr)) +
  # Shaded region showing significant interaction effect
  geom_rect(
    inherit.aes = FALSE,
    xmin = jn_lower, xmax = jn_upper,
    ymin = -Inf, ymax = Inf,
    fill = "gray80", alpha = 0.4
  ) +
  # Add actual RT lines
  geom_line(linewidth = 1.2) +
  # Add data points
  geom_point(size = 2, alpha = 0.8) +
  # Labels and styling
  labs(
    title = "Observed RT across Optimism and Accuracy",
    subtitle = paste0("Shaded region = Johnson-Neyman interval [", round(jn_lower, 2), ", ", round(jn_upper, 2), "]"),
    x = "Optimism (centered)",
    y = "Observed RT (ms)",
    colour = "Accuracy",
    linetype = "Accuracy"
  ) +
  scale_color_manual(values = c("red", "darkcyan")) +
  scale_linetype_manual(values = c("solid", "dashed")) +
  theme_minimal(base_size = 14)

## Warning: Removed 14 rows containing missing values or values outside the scale range
## (`geom_rect()`).

# Results Behaviour on Trial N+1

Post Response RT

Descriptives

#  Align ID types before joining
data.main.postbeh$ID <- as.character(data.main.postbeh$ID)
data.pers.cen$ID <- as.character(data.pers.cen$ID)

#  Calculate RT on trial N+1
data.main.postbeh <- data.main.postbeh %>%
  arrange(ID, block, trial_index) %>%
  group_by(ID) %>%
  mutate(RT_trial_nplus = lead(RT)) %>%
  ungroup()

#  Remove existing opt_cen or worry_cen columns (with extensions)
data.main.postbeh <- data.main.postbeh %>%
  select(-matches("^opt_cen"), -matches("^worry_cen"))

data.main.postbeh <- data.main.postbeh %>%
  left_join(
    data.pers.cen %>% select(ID, opt_cen, worry_cen),
    by = "ID"
  )
#  Compute mean RT on trial N+1 grouped by previous accuracy
matrix.mean.PES <- data.main.postbeh %>%
  filter(RT_trial_nplus != 999999) %>%
  group_by(ID, accuStr) %>%
  summarise(
    m.RTnplus = mean(RT_trial_nplus, na.rm = TRUE),
    se.RTnplus = sd(RT_trial_nplus, na.rm = TRUE) / sqrt(n()),
    opt_cen = dplyr::first(opt_cen),
    worry_cen = dplyr::first(worry_cen),
    .groups = "drop"
  )

#  Keep only participants who have both error & correct trials
valid_ids <- matrix.mean.PES %>%
  group_by(ID) %>%
  summarise(n = n()) %>%
  filter(n == 2) %>%
  pull(ID)

matrix.mean.PES <- matrix.mean.PES %>%
  filter(ID %in% valid_ids)

# factor order
matrix.mean.PES$accuStr <- factor(matrix.mean.PES$accuStr, levels = c("error", "correct"))

# Step 7: Pivot to wide format and compute PES
PES_data <- matrix.mean.PES %>%
  select(ID, accuStr, m.RTnplus, opt_cen, worry_cen) %>%
  pivot_wider(names_from = accuStr, values_from = m.RTnplus) %>%
  mutate(
    PES = error - correct  # Post-error slowing
  )

# Model PES ~ optimism and worry

summary(lm(PES ~ opt_cen + worry_cen, data = PES_data))

## 
## Call:
## lm(formula = PES ~ opt_cen + worry_cen, data = PES_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -206.393  -30.560   -1.623   34.235  211.620 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)    4.116      6.783   0.607   0.5455  
## opt_cen        5.456     12.826   0.425   0.6716  
## worry_cen    -12.910      7.710  -1.675   0.0976 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 64.71 on 88 degrees of freedom
## Multiple R-squared:  0.03118,    Adjusted R-squared:  0.009161 
## F-statistic: 1.416 on 2 and 88 DF,  p-value: 0.2481

# Summary of RT by accuracy
average.PES.Acc <- matrix.mean.PES %>%
  group_by(accuStr) %>%
  summarise(
    m.PES = mean(m.RTnplus, na.rm = TRUE),
    se.PES = sd(m.RTnplus, na.rm = TRUE) / sqrt(n()),
    .groups = "drop"
  )

# Summary of PES overall
average.PES.EV <- PES_data %>%
  summarise(
    meanRT = mean(PES, na.rm = TRUE),
    seRT = sd(PES, na.rm = TRUE) / sqrt(n())
  )

#  Print summary
cat(paste0("Post-error slowing:\n",
       "Error RT = ", round(average.PES.Acc$m.PES[average.PES.Acc$accuStr == "error"], 0), " ± ", round(average.PES.Acc$se.PES[average.PES.Acc$accuStr == "error"], 2), " ms\n",
       "Correct RT = ", round(average.PES.Acc$m.PES[average.PES.Acc$accuStr == "correct"], 0), " ± ", round(average.PES.Acc$se.PES[average.PES.Acc$accuStr == "correct"], 2), " ms\n",
       "PES (error – correct) = ", round(average.PES.EV$meanRT, 0), " ± ", round(average.PES.EV$seRT, 2), " ms\n"))

## Post-error slowing:
## Error RT = 894 ± 10.47 ms
## Correct RT = 890 ± 8.63 ms
## PES (error – correct) = 4 ± 6.81 ms

# Visual Diagnostics ===

# 1. Violin Plot (RT on trial n+1 by previous accuracy)
PES_long <- matrix.mean.PES %>%
  mutate(PES = m.RTnplus)

ggplot(PES_long, aes(x = accuStr, y = PES, color = accuStr)) +
  geom_violin(position = position_dodge(1), fill = NA) +
  geom_boxplot(position = position_dodge(1), width = 0.1, outlier.shape = NA) +
  geom_point(position = position_jitterdodge(jitter.width = 0.1, dodge.width = 1)) +

  # ✅ Add mean ± SE (confidence bands)
  stat_summary(fun = mean, geom = "point", shape = 23, size = 3, fill = "white", color = "black") +
  stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.1, color = "black") +

  scale_color_manual(values = c("#f20a19", "#3df50a")) +
  labs(
    x = "Accuracy Type (accuStr)",
    y = "RT on Trial N+1",
    color = "Accuracy"
  ) +
  theme_classic()

# 2. Histogram of participant-level PES
ggplot(PES_data, aes(x = PES)) +
  geom_histogram(bins = 30, fill = "steelblue", color = "white") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "red") +
  labs(title = "Distribution of Post-Error Slowing (PES) across Participants",
       x = "PES (RT after error – RT after correct)",
       y = "Count") +
  theme_minimal()

# 3. RT on Trial N+1 by Previous Accuracy (full distribution)
ggplot(data.main.postbeh, aes(x = RT_trial_nplus, fill = accuStr)) +
  geom_histogram(position = "identity", alpha = 0.5, bins = 50) +
  facet_wrap(~accuStr) +
  labs(title = "RT on Trial N+1 by Accuracy of Previous Trial",
       x = "RT (ms)", y = "Count") +
  theme_minimal()

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

LME - PES, the famous 4way interaction?

##1. Visualize RT on trial N+1 grouped by previous trial accuracy (accuStr)

# RT on trial N+1
data.main <- data.main %>%
  arrange(ID, trial_index) %>%
  group_by(ID) %>%
  mutate(RT_trial_nplus = lead(RT)) %>%
  ungroup()

# Calculate mean RT_n+1 for each participant × accuStr
matrix.mean.PES_long <- data.main %>%
  filter(!is.na(RT_trial_nplus)) %>%
  group_by(ID, accuStr) %>%
  summarise(
    PES = mean(RT_trial_nplus, na.rm = TRUE),
    opt_cen = first(opt_cen),
    worry_cen = first(worry_cen),
    .groups = "drop"
  )

# Plot violin for trial N+1 RT by accuracy on trial N
ggplot(matrix.mean.PES_long, aes(x = accuStr, y = PES, color = accuStr)) +
  geom_violin(fill = NA) +
  geom_boxplot(width = 0.1, outlier.shape = NA) +
  geom_point(position = position_jitter(width = 0.1), alpha = 0.7) +
  labs(x = "Previous Trial Accuracy (accuStr)", y = "RT on Trial N+1 [ms]") +
  theme_classic()

##2. Compute a single PES score per participant for modeling

#For regression: lm(PES ~ opt_cen * worry_cen)
# Create label for trial type (N)
data.main <- data.main %>%
  mutate(Type_Nplus = case_when(
    accuStr == "error" ~ "post_error",
    accuStr == "correct" ~ "post_correct",
    TRUE ~ NA_character_
  ))

# Compute PES as RT_post-error – RT_post-correct
matrix.mean.PES <- data.main %>%
  filter(!is.na(RT_trial_nplus)) %>%
  group_by(ID, Type_Nplus) %>%
  summarise(
    mean_RT = mean(RT_trial_nplus, na.rm = TRUE),
    opt_cen = first(opt_cen),
    worry_cen = first(worry_cen),
    .groups = "drop"
  ) %>%
  pivot_wider(names_from = Type_Nplus, values_from = mean_RT) %>%
  mutate(PES = post_error - post_correct)

# Fit regression model
model_pes <- lm(PES ~ opt_cen * worry_cen, data = matrix.mean.PES)
summary(model_pes)

## 
## Call:
## lm(formula = PES ~ opt_cen * worry_cen, data = matrix.mean.PES)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -207.969  -30.655   -2.613   33.282  212.719 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)
## (Intercept)          4.659      7.213   0.646    0.520
## opt_cen              5.273     12.920   0.408    0.684
## worry_cen          -12.854      7.756  -1.657    0.101
## opt_cen:worry_cen   -2.947     12.750  -0.231    0.818
## 
## Residual standard error: 65.06 on 87 degrees of freedom
## Multiple R-squared:  0.03177,    Adjusted R-squared:  -0.001613 
## F-statistic: 0.9517 on 3 and 87 DF,  p-value: 0.4194

# Optional plots
ggplot(matrix.mean.PES, aes(x = opt_cen, y = PES)) +
  geom_point() +
  geom_smooth(method = "lm") +
  labs(title = "PES as a Function of Optimism",
       x = "Optimism (Centered)", y = "Post-Error Slowing (ms)") +
  theme_minimal()

## `geom_smooth()` using formula = 'y ~ x'

ggplot(matrix.mean.PES, aes(x = worry_cen, y = PES)) +
  geom_point() +
  geom_smooth(method = "lm") +
  labs(title = "PES as a Function of Worry",
       x = "Worry (Centered)", y = "Post-Error Slowing (ms)") +
  theme_minimal()

## `geom_smooth()` using formula = 'y ~ x'

ggplot(matrix.mean.PES_long, aes(x = accuStr, y = PES, color = accuStr)) +
  geom_violin(fill = NA) +
  geom_boxplot(width = 0.1, outlier.shape = NA) +
  geom_point(position = position_jitter(width = 0.1), alpha = 0.7) +
  labs(
    x = "Previous Trial Accuracy (accuStr)",
    y = "RT on Trial N+1 [ms]",
    color = "accuStr"
  ) +
  theme_classic()

Collinearity - PES

check_collinearity(model_pes)

## # Check for Multicollinearity
## 
## Low Correlation
## 
##               Term  VIF       VIF 95% CI adj. VIF Tolerance Tolerance 95% CI
##            opt_cen 1.14 [1.03,     1.80]     1.07      0.88     [0.56, 0.98]
##          worry_cen 1.14 [1.02,     1.81]     1.07      0.88     [0.55, 0.98]
##  opt_cen:worry_cen 1.00 [1.00,      Inf]     1.00      1.00     [0.00, 1.00]

collinearity_result <- check_collinearity(model_pes)
print(collinearity_result)

## # Check for Multicollinearity
## 
## Low Correlation
## 
##               Term  VIF       VIF 95% CI adj. VIF Tolerance Tolerance 95% CI
##            opt_cen 1.14 [1.03,     1.80]     1.07      0.88     [0.56, 0.98]
##          worry_cen 1.14 [1.02,     1.81]     1.07      0.88     [0.55, 0.98]
##  opt_cen:worry_cen 1.00 [1.00,      Inf]     1.00      1.00     [0.00, 1.00]

# plot results
if (require("see")) {
  x <- check_collinearity(model_pes)
  plot(x)
}

Post hoc Acc x Eval type

model_pes_full <- lm(PES ~ opt_cen * worry_cen, data = matrix.mean.PES)
summary(model_pes_full)

## 
## Call:
## lm(formula = PES ~ opt_cen * worry_cen, data = matrix.mean.PES)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -207.969  -30.655   -2.613   33.282  212.719 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)
## (Intercept)          4.659      7.213   0.646    0.520
## opt_cen              5.273     12.920   0.408    0.684
## worry_cen          -12.854      7.756  -1.657    0.101
## opt_cen:worry_cen   -2.947     12.750  -0.231    0.818
## 
## Residual standard error: 65.06 on 87 degrees of freedom
## Multiple R-squared:  0.03177,    Adjusted R-squared:  -0.001613 
## F-statistic: 0.9517 on 3 and 87 DF,  p-value: 0.4194

ggplot(matrix.mean.PES, aes(x = opt_cen, y = PES, color = worry_cen)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  theme_minimal()

## `geom_smooth()` using formula = 'y ~ x'

## Warning: The following aesthetics were dropped during statistical transformation:
## colour.
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
##   the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
##   variable into a factor?

colnames(PES_data)

## [1] "ID"        "opt_cen"   "worry_cen" "correct"   "error"     "PES"

Figures (Appendix)

matrix.mean.PES <- matrix.mean.PES %>%
  mutate(accuStr = "delta")

ggplot(matrix.mean.PES, aes(x = "delta", y = PES)) +
  geom_violin(fill = "skyblue", alpha = 0.3) +
  geom_boxplot(width = 0.1, fill = "white") +
  geom_point(position = position_jitter(width = 0.1), alpha = 0.7) +
  stat_summary(fun = mean, geom = "crossbar", width = 0.3, color = "black") +
  labs(x = "", y = "Post-Error Slowing (PES) [ms]") +
  theme_classic()

Post Response accuracy {.tabset .tabset-pills}, triple check*

Descriptives

# prepare table
matrix.mean.PEA <- data.main.postbeh  %>%
  group_by(ID, accuStr ) %>%
  summarise(correct = mean(correct_nplus)*100,
            opt_cen = mean(opt_cen),
            worry_cen = mean(worry_cen)
            ) %>%
  ungroup()

# Difference 
matrix.delta.PEA <- matrix.mean.PEA %>%
   group_by(ID) %>%
   summarise(correct = diff(correct),
            opt_cen = first(opt_cen),
            worry_cen = first(worry_cen)) %>%
  ungroup()  

data.main.postbeh <- data.main.postbeh %>%
  mutate(
    accuNum = ifelse(accuStr == "correct", -0.5, 0.5),
    
  )
matrix.mean.PEA.b <- data.main.postbeh %>%
  group_by(ID, accuNum) %>%
  summarise(
    correct = mean(correct_nplus, na.rm = TRUE) * 100,
    opt_cen = mean(opt_cen, na.rm = TRUE),
    worry_cen = mean(worry_cen, na.rm = TRUE)
  ) %>%
  ungroup()


PEA.c.lmeb <- lmer(
  correct ~ accuNum * opt_cen * worry_cen + (1 | ID),
  data = matrix.mean.PEA.b
)

# Define standard error function
stderror <- function(x) {
  x <- x[!is.na(x)]
  sd(x) / sqrt(length(x))
}

# Compute mean and SE by accuracy type
average.PEA.Acc <- matrix.mean.PEA %>%
  group_by(accuStr) %>%
  summarise(
    m.correct = mean(correct, na.rm = TRUE),
    se.correct = stderror(correct),
    .groups = "drop"
  )

tab_model(PEA.c.lmeb, dv.labels = "Post-Response Accuracy")

	Post-Response Accuracy
Predictors	Estimates	CI	p
(Intercept)	84.76	82.41 – 87.11	<0.001
accuNum	0.11	-1.66 – 1.88	0.900
opt cen	0.20	-4.01 – 4.41	0.926
worry cen	-0.85	-3.38 – 1.67	0.506
accuNum × opt cen	-1.38	-4.55 – 1.79	0.392
accuNum × worry cen	1.27	-0.64 – 3.17	0.191
opt cen × worry cen	-0.63	-4.79 – 3.52	0.765
(accuNum × opt cen) × worry cen	-0.54	-3.67 – 2.59	0.736
Random Effects
σ2	32.73
τ00 ID	99.04
ICC	0.75
N ID	91
Observations	182
Marginal R2 / Conditional R2	0.008 / 0.754

paste0("Post response accuracy: ",
       paste(average.PEA.Acc$accuStr, 
             round(average.PEA.Acc$m.correct, 2), 
             " ± ", 
             round(average.PEA.Acc$se.correct, 2), 
             "%", 
             collapse = "; "))

## [1] "Post response accuracy: correct 78.48  ±  3.9 %; error 84.71  ±  1.43 %"

Statistical analyses

LME - PEA

PEA.c.lm <- lm(correct ~ accuStr * opt_cen * worry_cen, data = matrix.mean.PEA)
summary(PEA.c.lm)

## 
## Call:
## lm(formula = correct ~ accuStr * opt_cen * worry_cen, data = matrix.mean.PEA)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -44.971  -7.044   3.175   9.335  16.353 
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      78.969      3.810  20.726   <2e-16 ***
## accuStrerror                      5.932      4.116   1.441    0.153    
## opt_cen                         -10.350     15.466  -0.669    0.505    
## worry_cen                        -4.944      3.502  -1.412    0.161    
## accuStrerror:opt_cen              9.967     15.706   0.635    0.527    
## accuStrerror:worry_cen            4.697      3.883   1.209    0.230    
## opt_cen:worry_cen               -15.270     15.126  -1.010    0.315    
## accuStrerror:opt_cen:worry_cen   14.287     15.364   0.930    0.355    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13.65 on 92 degrees of freedom
##   (82 observations deleted due to missingness)
## Multiple R-squared:  0.0571, Adjusted R-squared:  -0.01465 
## F-statistic: 0.7958 on 7 and 92 DF,  p-value: 0.5927

# binary - required for JN-procedure
PEA.c.lmeb <- lmer(correct ~ accuNum * opt_cen * worry_cen + (1  | ID), matrix.mean.PEA.b)

lm(correct ~ accuStr * opt_cen + accuStr * worry_cen + opt_cen * worry_cen, data = matrix.mean.PEA)

## 
## Call:
## lm(formula = correct ~ accuStr * opt_cen + accuStr * worry_cen + 
##     opt_cen * worry_cen, data = matrix.mean.PEA)
## 
## Coefficients:
##            (Intercept)            accuStrerror                 opt_cen  
##                 78.809                   6.178                  -4.827  
##              worry_cen    accuStrerror:opt_cen  accuStrerror:worry_cen  
##                 -4.770                   4.416                   4.534  
##      opt_cen:worry_cen  
##                 -1.421

lm(correct ~ accuStr + opt_cen + worry_cen, data = matrix.mean.PEA)

## 
## Call:
## lm(formula = correct ~ accuStr + opt_cen + worry_cen, data = matrix.mean.PEA)
## 
## Coefficients:
##  (Intercept)  accuStrerror       opt_cen     worry_cen  
##    78.564101      6.139520     -0.002958     -1.091992

#Visualizing Interactions
interact_plot(
  model = lm(correct ~ accuStr + opt_cen + worry_cen, data = matrix.mean.PEA),
  pred = "opt_cen",
  modx = "accuStr",
  plot.points = TRUE
)

## Warning: opt_cen and accuStr are not included in an interaction with one another in
## the model.

#Scatter-plot
ggplot(matrix.mean.PEA, aes(x = opt_cen, y = correct)) +
  geom_point() +
  geom_smooth(method = "lm") +
  labs(
    title = "Accuracy vs. Optimism",
    x = "Optimism (centered)",
    y = "Post-response Accuracy (%)"
  ) +
  theme_minimal()

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 82 rows containing non-finite outside the scale range
## (`stat_smooth()`).

## Warning: Removed 82 rows containing missing values or values outside the scale range
## (`geom_point()`).

#JN Plot
johnson_neyman(
  model = PEA.c.lmeb,
  pred = "opt_cen",
  modx = "accuNum",
  mod.range = c(-0.5, 0.5)
)

## JOHNSON-NEYMAN INTERVAL
## 
## The Johnson-Neyman interval could not be found. Is the p value for your
## interaction term below the specified alpha?

# Assuming matrix.mean.PEA.b contains 'opt_cen' and 'correct'
ggplot(matrix.mean.PEA.b, aes(x = opt_cen, y = correct)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "lm", color = "steelblue", fill = "lightgray") +
  labs(
    title = "Effect of Optimism on Post-Response Accuracy",
    x = "Optimism (centered)",
    y = "Post-response Accuracy (%)"
  ) +
  theme_classic()

## `geom_smooth()` using formula = 'y ~ x'

Collinearity - PEA

check_collinearity(PEA.c.lm)

## Model has interaction terms. VIFs might be inflated.
##   Try to center the variables used for the interaction, or check
##   multicollinearity among predictors of a model without interaction terms.

## # Check for Multicollinearity
## 
## Low Correlation
## 
##     Term  VIF     VIF 95% CI adj. VIF Tolerance Tolerance 95% CI
##  accuStr 1.03 [ 1.00, 24.40]     1.01      0.97     [0.04, 1.00]
## 
## Moderate Correlation
## 
##               Term  VIF     VIF 95% CI adj. VIF Tolerance Tolerance 95% CI
##          worry_cen 6.03 [ 4.42,  8.40]     2.46      0.17     [0.12, 0.23]
##  accuStr:worry_cen 6.17 [ 4.52,  8.60]     2.48      0.16     [0.12, 0.22]
## 
## High Correlation
## 
##                       Term   VIF     VIF 95% CI adj. VIF Tolerance
##                    opt_cen 38.09 [26.98, 53.96]     6.17      0.03
##            accuStr:opt_cen 38.06 [26.95, 53.92]     6.17      0.03
##          opt_cen:worry_cen 33.34 [23.63, 47.21]     5.77      0.03
##  accuStr:opt_cen:worry_cen 33.26 [23.58, 47.10]     5.77      0.03
##  Tolerance 95% CI
##      [0.02, 0.04]
##      [0.02, 0.04]
##      [0.02, 0.04]
##      [0.02, 0.04]

# plot results
if (require("see")) {
  x <- check_collinearity(PEA.c.lm)
  plot(x)
}

## Model has interaction terms. VIFs might be inflated.
##   Try to center the variables used for the interaction, or check
##   multicollinearity among predictors of a model without interaction terms.

Post hoc simple slope - OPT x Accu

#Newest code 
# One-time reassignment of accuStr from accuNum
matrix.mean.PEA.b <- matrix.mean.PEA.b %>%
  mutate(accuStr = factor(ifelse(accuNum == -0.5, "correct", "error")))

# Clean simple slope plot
ggplot(matrix.mean.PEA.b, aes(x = opt_cen, y = correct, color = accuStr, linetype = accuStr)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "lm", se = TRUE, alpha = 0.1) +
  
  # Optional JN range shading (adjust bounds if needed)
  geom_rect(
    inherit.aes = FALSE,
    fill = "lightgray",
    alpha = 0.2,
    xmin = -0.2, xmax = 0.4,
    ymin = -Inf, ymax = Inf
  ) +
  
  labs(
    title = "Simple Slopes: Optimism Predicting Accuracy by Previous Trial",
    x = "Optimism (centered)",
    y = "Post-Response Accuracy (%)",
    color = "Previous Accuracy",
    linetype = "Previous Accuracy"
  ) +
  scale_color_grey(start = 0, end = 0.2) +
  scale_linetype_manual(values = c("solid", "longdash")) +
  theme_classic()

## `geom_smooth()` using formula = 'y ~ x'

# Run separate models for each condition
lm(correct ~ opt_cen, data = matrix.mean.PEA.b %>% filter(accuStr == "correct"))

## 
## Call:
## lm(formula = correct ~ opt_cen, data = matrix.mean.PEA.b %>% 
##     filter(accuStr == "correct"))
## 
## Coefficients:
## (Intercept)      opt_cen  
##    84.63992      0.04573

lm(correct ~ opt_cen, data = matrix.mean.PEA.b %>% filter(accuStr == "error"))

## 
## Call:
## lm(formula = correct ~ opt_cen, data = matrix.mean.PEA.b %>% 
##     filter(accuStr == "error"))
## 
## Coefficients:
## (Intercept)      opt_cen  
##     84.6546      -0.5729

###########

# Archived Code
# 1. Compute prediction grid
limits <- max(abs(c(
  range(data.main$opt_cen, na.rm = TRUE),
  range(data.main$worry_cen, na.rm = TRUE)
)))

b <- seq(-0.5 - limits, 0.5 + limits, by = 0.1)

data.pred <- expand.grid(
  ID = unique(data.main$ID),
  accuStr = c("correct", "error"),
  worry_cen = 0,  # hold constant
  opt_cen = b
) %>%
  mutate(
    accuStr = factor(accuStr, levels = c("correct", "error")),
    accuNum = ifelse(accuStr == "correct", -0.5, 0.5)
  )

# 2. Predict using lmer model
data.pred$PEA_correct <- predict(PEA.c.lmeb, newdata = data.pred, re.form = NA)


# 4. Plot
matrix.mean.PEA.b <- matrix.mean.PEA.b %>%
  mutate(
    accuStr = factor(ifelse(accuNum == -0.5, "correct", "error"))
  )

library(ggplot2)

ggplot(matrix.mean.PEA.b, aes(x = opt_cen, y = correct, color = accuStr, linetype = accuStr)) +
  geom_smooth(method = "lm", se = TRUE, alpha = 0.1) +
  labs(
    x = "Optimism (centered)",
    y = "Post-Response Accuracy (%)",
    color = "Previous Accuracy",
    linetype = "Previous Accuracy"
  ) +
  scale_color_grey(start = 0, end = 0.2) +
  scale_linetype_manual(values = c("solid", "longdash")) +
  theme_classic()

## `geom_smooth()` using formula = 'y ~ x'

ggplot(matrix.mean.PEA.b, aes(x = opt_cen, y = correct, color = as.factor(accuNum))) +
  geom_point() +
  geom_smooth(method = "lm") +
  labs(
    title = "Post-response Accuracy by Optimism and Previous Accuracy",
    x = "Optimism (centered)",
    y = "Post-response Accuracy (%)",
    color = "Previous Accuracy"
  ) +
  theme_minimal()

## `geom_smooth()` using formula = 'y ~ x'

Post hoc simple slope - OPT x Accu (again)

# Clean version (run mutate only once)
matrix.mean.PEA.b <- matrix.mean.PEA.b %>%
  mutate(accuStr = factor(ifelse(accuNum == -0.5, "correct", "error")))

# Plot
ggplot(matrix.mean.PEA.b, aes(x = opt_cen, y = correct, color = accuStr, linetype = accuStr)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "lm", se = TRUE, alpha = 0.1) +
  geom_rect(
    inherit.aes = FALSE,
    fill = "lightgray",
    alpha = 0.2,
    xmin = -0.2, xmax = 0.4,  # Adjust this manually based on JN output
    ymin = -Inf, ymax = Inf
  ) +
  labs(
    title = "Simple Slopes: Optimism Predicting Post-Response Accuracy",
    x = "Optimism (centered)",
    y = "Accuracy (%)",
    color = "Previous Accuracy",
    linetype = "Previous Accuracy"
  ) +
  theme_classic()

## `geom_smooth()` using formula = 'y ~ x'

geom_rect(
  inherit.aes = FALSE,
  fill = "lightgray",
  alpha = 0.2,
  xmin = -0.2, xmax = 0.4,
  ymin = -Inf, ymax = Inf
)

## geom_rect: linejoin = mitre, na.rm = FALSE
## stat_identity: na.rm = FALSE
## position_identity

# Report simple slopes
lm(correct ~ opt_cen, data = matrix.mean.PEA.b %>% filter(accuStr == "correct"))

## 
## Call:
## lm(formula = correct ~ opt_cen, data = matrix.mean.PEA.b %>% 
##     filter(accuStr == "correct"))
## 
## Coefficients:
## (Intercept)      opt_cen  
##    84.63992      0.04573

lm(correct ~ opt_cen, data = matrix.mean.PEA.b %>% filter(accuStr == "error"))

## 
## Call:
## lm(formula = correct ~ opt_cen, data = matrix.mean.PEA.b %>% 
##     filter(accuStr == "error"))
## 
## Coefficients:
## (Intercept)      opt_cen  
##     84.6546      -0.5729

Figures (Appendix)

# Fig. RT interaction -- mean distributions  
  col4 <- c("#3df50a", "#f20a19")  # green and red for correct/error

ggplot(matrix.mean.RT, aes(x = accuStr, y = mRT, color = accuStr)) +
  geom_violin(position = position_dodge(1), fill = NA) +
  geom_boxplot(position = position_dodge(1), width = 0.1, outlier.shape = NA) +
  geom_point(position = position_jitterdodge(jitter.width = 0.1, dodge.width = 1), alpha = 0.7) +
  scale_color_manual(values = col4) +
  labs(
    x = "Previous Trial Accuracy",
    y = "Response Time [ms]",
    color = "Accuracy"
  ) +
  theme_classic()

##block_wise

data.correct <- data.main.postbeh %>%
  filter(accuStr == "correct")

matrix.RT.block.correct <- data.correct %>%
  group_by(ID, block) %>%
  summarise(
    m.RT = mean(RT, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  mutate(blockNEW = as.factor(block))  
ggplot(matrix.RT.block.correct, aes(x = blockNEW, y = m.RT)) +
  geom_violin(fill = "lightblue", alpha = 0.3) +
  geom_boxplot(width = 0.1, outlier.shape = NA) +
  geom_point(position = position_jitter(width = 0.1), alpha = 0.6) +
  labs(
    x = "Block",
    y = "Response Time [ms]"
  ) +
  theme_classic()

Results Elecrotrophysiology

Error negativity (Ne /c)

Descriptives

stderror <- function(x) {
  x <- x[!is.na(x)]
  sd(x) / sqrt(length(x))
}

# Grand average: Ne by accuracy
average.Ne4.Acc <- matrix.mean.erp %>%
  group_by(acc) %>%
  summarise(
    m.Ne = mean(Ne, na.rm = TRUE),
    se.Ne = stderror(Ne)
  ) %>%
  ungroup()

# Grand average: overall Ne
average.Ne4.oa <- matrix.mean.erp %>%
  summarise(
    m.Ne = mean(Ne, na.rm = TRUE),
    se.Ne = stderror(Ne)
  )

# Print result
paste("Mean Ne by accuracy:",
      paste(average.Ne4.Acc$acc, round(average.Ne4.Acc$m.Ne, 3), "+/-", round(average.Ne4.Acc$se.Ne, 3), "µV/cm²", collapse = "; "))

## [1] "Mean Ne by accuracy: correct -14.855 +/- 1.002 µV/cm²; error -24.793 +/- 1.793 µV/cm²"

Statistical analyses

LME - NE***triple check

Ne4.0.lme <- lmer(Ne ~ acc *  opt_cen  * worry_cen + (1 | ID), matrix.mean.erp)
tab_summary(Ne4.0.lme, pred.labels = predictors.all, dv.labels = "Ne/c Amplitude Effects")

	Ne/c Amplitude Effects
Predictors	Estimate	SE	df	t	p
(Intercept)	-19.56	1.29	87.00	-15.11	<0.001
acc1	-9.43	1.65	87.00	-5.72	<0.001
Optimism	-3.16	2.32	87.00	-1.36	0.177
Worry	2.18	1.39	87.00	1.57	0.120
acc1:opt_cen	-2.01	2.95	87.00	-0.68	0.499
acc1:worry_cen	2.23	1.77	87.00	1.26	0.212
Response Type x Worry	-1.44	2.29	87.00	-0.63	0.529
acc1:opt_cen:worry_cen	-2.77	2.91	87.00	-0.95	0.345
Random Effects
σ2	110.50
τ00 ID	81.03
ICC	0.42
N ID	91
Observations	182
Marginal R2 / Conditional R2	0.142 / 0.505

Collinearity

check_collinearity(Ne4.0.lme )

## # Check for Multicollinearity
## 
## Low Correlation
## 
##                   Term  VIF       VIF 95% CI adj. VIF Tolerance
##                    acc 1.12 [1.03,     1.50]     1.06      0.89
##                opt_cen 1.14 [1.04,     1.49]     1.07      0.88
##              worry_cen 1.14 [1.04,     1.49]     1.07      0.88
##            acc:opt_cen 1.14 [1.04,     1.49]     1.07      0.88
##          acc:worry_cen 1.14 [1.04,     1.49]     1.07      0.88
##      opt_cen:worry_cen 1.00 [1.00, 9.71e+12]     1.00      1.00
##  acc:opt_cen:worry_cen 1.12 [1.03,     1.49]     1.06      0.89
##  Tolerance 95% CI
##      [0.67, 0.97]
##      [0.67, 0.96]
##      [0.67, 0.96]
##      [0.67, 0.96]
##      [0.67, 0.96]
##      [0.00, 1.00]
##      [0.67, 0.97]

# plot results
if (require("see")) {
  x <- check_collinearity(Ne4.0.lme)
  plot(x)
}

Post hoc simple slope Acc x OPT x WORRY

# --- POST HOC: 3-way Interaction Plot (Optimism × Worry × Accuracy Effect on Ne) ---
library(lme4)

Ne3.lme <- lmer(Ne ~ acc * opt_cen * worry_cen + (1 | ID), data = matrix.mean.erp)

# Generate interaction surface using lab-specific function (you should already have it)
ssp3 <- list(
  Ne3 = sim_slope_plot_3x(Ne3.lme, steps = 0.05)
)

# Customize plot appearance
ssp3$Ne3 <- ssp3$Ne3 +
  scale_fill_gradient2(
    limits = c(-0.3, 0.3),
    low = muted("blue"),
    mid = "white",
    high = muted("red"),
    name = "Accuracy Effect"
  ) +
  theme(legend.position = "bottom") +
  labs(x = "Optimism (centered)", y = "Worry (centered)")

## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.

# Show the final plot
cowplot::plot_grid(ssp3$Ne3)

#Optimism × Accuracy interaction only (Ne amplitude)
Ne_opt.lme <- lmer(m.Ne ~ acc * opt_cen + (1 | ID), data = matrix.mean.Ne3)

interact_plot(Ne3.lme, pred = "opt_cen", modx = "acc", 
              interval = TRUE, plot.points = TRUE,
              y.label = "Ne Amplitude (µV/cm²)",
              x.label = "Optimism (centered)",
              modx.labels = c("Correct", "Error")) +
  theme_classic()

Error positivity (Pe)

Descriptives

matrix.mean.Pe2 <- matrix.mean.erp %>%
  group_by(ID) %>%
  summarise(m.Pe = mean(Pe),
            se.Pe = sd(Pe),
            opt_cen =first(opt_cen),
            worry_cen = first(worry_cen)) %>%
  ungroup()


#grand average :  accuracy 
 average.Pe.Acc <- matrix.mean.erp  %>%
  group_by(acc) %>%
  summarise(m.Pe = mean(Pe),
            se.Pe = stderror(Pe)) %>%
  ungroup()   
  

paste("mean Pe/c:", average.Pe.Acc$acc, round(average.Pe.Acc$m.Pe, 3), "+/-",round(average.Pe.Acc$se.Pe, 3), "µV/cm²")

## [1] "mean Pe/c: correct 14.231 +/- 1.57 µV/cm²"
## [2] "mean Pe/c: error 32.827 +/- 2.087 µV/cm²"

Statistical Analyses

LME - Pe

Pe4.lme <- lmer(Pe ~ acc *  opt_cen  * worry_cen + (1 | ID), matrix.mean.erp)
tab_summary(Pe4.lme, pred.labels = predictors.all, dv.labels = "Pe/c Amplitude Effects")

	Pe/c Amplitude Effects
Predictors	Estimate	SE	df	t	p
(Intercept)	23.44	1.65	87.00	14.17	<0.001
acc1	19.13	2.14	87.00	8.96	<0.001
Optimism	-0.88	2.96	87.00	-0.30	0.766
Worry	-1.75	1.78	87.00	-0.99	0.327
acc1:opt_cen	2.21	3.83	87.00	0.58	0.566
acc1:worry_cen	-2.50	2.30	87.00	-1.09	0.280
Response Type x Worry	0.49	2.92	87.00	0.17	0.867
acc1:opt_cen:worry_cen	-2.91	3.78	87.00	-0.77	0.444
Random Effects
σ2	185.65
τ00 ID	129.63
ICC	0.41
N ID	91
Observations	182
Marginal R2 / Conditional R2	0.227 / 0.545

Collinearity

check_collinearity(Pe4.lme)

## # Check for Multicollinearity
## 
## Low Correlation
## 
##                   Term  VIF       VIF 95% CI adj. VIF Tolerance
##                    acc 1.12 [1.03,     1.50]     1.06      0.89
##                opt_cen 1.14 [1.04,     1.49]     1.07      0.88
##              worry_cen 1.14 [1.04,     1.49]     1.07      0.88
##            acc:opt_cen 1.14 [1.04,     1.49]     1.07      0.88
##          acc:worry_cen 1.14 [1.04,     1.49]     1.07      0.88
##      opt_cen:worry_cen 1.00 [1.00, 9.71e+12]     1.00      1.00
##  acc:opt_cen:worry_cen 1.12 [1.03,     1.49]     1.06      0.89
##  Tolerance 95% CI
##      [0.67, 0.97]
##      [0.67, 0.96]
##      [0.67, 0.96]
##      [0.67, 0.96]
##      [0.67, 0.96]
##      [0.00, 1.00]
##      [0.67, 0.97]

# plot results
if (require("see")) {
  x <- check_collinearity(Pe4.lme)
  plot(x)
}

Posthoc Acc

# within interaction post hoc comparison
  emmeans(Pe4.lme, specs = c("acc")) %>% 
  pairs()

## NOTE: Results may be misleading due to involvement in interactions

##  contrast        estimate   SE df t.ratio p.value
##  correct - error    -19.1 2.14 87  -8.955  <.0001
## 
## Degrees-of-freedom method: kenward-roger

Post hoc simple slope - Pe

Pe2.lm <- lm(m.Pe ~ opt_cen * worry_cen + (1 | ID), matrix.mean.Pe2)
tab_summary(Pe2.lm, pred.labels =  predictors.opt, dv.labels = "Pe/c Amplitude Effects")

## Model matrix is rank deficient. Parameters `1 | IDTRUE` were not
##   estimable.

## `ci_method` must be one of "residual", "wald", "normal", "profile",
##   "boot", "uniroot", "betwithin" or "ml1". Using "wald" now.

	Pe/c Amplitude Effects
Predictors	Estimate	SE	df	t	p
(Intercept)	23.44	1.65	87.00	14.17	<0.001
Optimism	-0.88	2.96	87.00	-0.30	0.766
Worry	-1.75	1.78	87.00	-0.99	0.327
Worry x Optimism	0.49	2.92	87.00	0.17	0.867
Observations	91
R2 / R2 adjusted	0.017 / -0.017

sspe2 <- list(Pe2= sim_slope_plot_2x(Pe2.lm, steps = 0.05, thresholding = TRUE, center = TRUE)) 

sspe2$Pe2 <- sspe2$Pe2 +
  scale_fill_gradient2(limits = c(-0.3,0.3), low = muted("blue"), mid = "white", high = muted("red"),  name = "Deviation from Mean Pe-Amplitue") +
  theme(legend.position = "bottom") + labs(x = "Optimism (centered)", y = "Worry (centered)")

## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.

cowplot::plot_grid(sspe2$Pe2 )

# Required libraries
library(ggplot2)

# Plot Pe amplitude by optimism (centered)
ggplot(matrix.mean.Pe2, aes(x = opt_cen, y = m.Pe)) +
  geom_point(color = "black", alpha = 0.6) +
  geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "blue") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "gray") +
  labs(
    title = "Pe Amplitude by Optimism (Centered)",
    x = "Optimism (centered)",
    y = "Pe Amplitude (µV/cm²)"
  ) +
  theme_classic()