WPT.Analyses.08.06.2020

Load packages

library(BBmisc)
library(corrr)
library(ggplot2)
library(lattice)
library(Hmisc)
library(corrplot)
library(tidyverse)
library(lme4)
library(sjPlot)
library(psych)
library(mediation)
library(lavaan)
library(pbkrtest)
library(Rcpp)
library(dplyr)

Load Data

WPT1 <- read.csv("wpt_study1_6.19.2020.csv", header=T, stringsAsFactors = FALSE, na.strings=c("","NA"))
WPT2 <- read.csv("Pred2Full.08.06.2020.csv", header=T, stringsAsFactors = FALSE, na.strings=c("","NA"))
WPT3 <- read.csv("Pred3Full.08.06.2020.csv", header=T, stringsAsFactors = FALSE, na.strings=c("","NA"))

show # of participants with below 50% accuracy and remove

#compute average accuracy for study 1
WPT1$acc <- as.integer(as.character(WPT1$acc))
averga_acc <- WPT1 %>%
  group_by(Participant)%>%
  dplyr::summarise(avg_acc=mean(acc, na.rm = T))

`summarise()` ungrouping output (override with `.groups` argument)

WPT1 <- merge(WPT1, averga_acc, by = "Participant", all.x = F)

#plot # of participants under 50% accurate
count_NAStudy1 <- WPT1$Participant[which(WPT1$avg_acc<=.52)]
count_NAStudy1<- unique(count_NAStudy1)
length(count_NAStudy1)

[1] 9

WPT1 <- WPT1[!(WPT1$Participant %in% count_NAStudy1),]

#plot # of participants under 50% accurate
count_NAStudy2 <- WPT2$Participant[which(WPT2$avg_acc<=.52)]
count_NAStudy2<- unique(count_NAStudy2)
length(count_NAStudy2)

[1] 20

WPT2 <- WPT2[!(WPT2$Participant %in% count_NAStudy2),] #remove those participants

#plot # of participants under 50% accurate
count_NAStudy3 <- WPT3$Participant[which(WPT3$avg_acc<=.52)]
count_NAStudy3<- unique(count_NAStudy3)
length(count_NAStudy3)

[1] 15

WPT3 <- WPT3[!(WPT3$Participant %in% count_NAStudy3),]

calculate average accuracy for each condition and plot

WPT1$acc <- as.numeric(as.character(WPT1$acc))
WPT1Summary <- WPT1 %>% # the names of the new data frame and the data frame to be summarised
  group_by(Condition) %>%   # the grouping variable
  dplyr::summarise(mean_acc = mean(acc, na.rm = T),  # calculates the mean of each group
            sd_PL = sd(acc),
            n_PL = n(),  # calculates the sample size per group
            SE_PL = sd(acc)/sqrt(n())) # calculates the standard error of each group

`summarise()` ungrouping output (override with `.groups` argument)

WPT2Summary <- WPT2 %>% # the names of the new data frame and the data frame to be summarised
  group_by(Condition) %>%   # the grouping variable
  dplyr::summarise(mean_acc = mean(acc, na.rm = T),  # calculates the mean of each group
            sd_PL = sd(acc),
            n_PL = n(),  # calculates the sample size per group
            SE_PL = sd(acc)/sqrt(n())) # calculates the standard error of each group

`summarise()` ungrouping output (override with `.groups` argument)

WPT3Summary <- WPT3 %>% # the names of the new data frame and the data frame to be summarised
  group_by(Condition) %>%# the grouping variable
  dplyr::summarise(meanAcc = mean(acc, na.rm = T),# calculates the mean of each group
                   sd_PL = sd(meanAcc),
                   n_PL = n(),# calculates the sample size per group
                   SE_PL = .01) # calculates the standard error of each group

`summarise()` ungrouping output (override with `.groups` argument)

#Study 1 
ggplot(WPT1Summary, aes(x=as.factor(Condition), y=mean_acc)) +
  geom_bar(stat = "identity", alpha=0.5) +
  geom_errorbar(aes(x=as.factor(Condition), ymin=mean_acc-SE_PL, ymax=mean_acc+SE_PL))+
  coord_cartesian(ylim = c(0, 1))

#Study 
ggplot(WPT2Summary, aes(x=as.factor(Condition), y=mean_acc)) +
  geom_bar(stat = "identity", alpha=0.5) +
  geom_errorbar(aes(x=as.factor(Condition), ymin=mean_acc-SE_PL, ymax=mean_acc+SE_PL))+
  coord_cartesian(ylim = c(0, 1))

#Study 2
ggplot(WPT3Summary, aes(x=as.factor(Condition), y=meanAcc)) +
  geom_bar(stat = "identity", alpha=0.5) +
  geom_errorbar(aes(x=as.factor(Condition), ymin=meanAcc-SE_PL, ymax=meanAcc+SE_PL))+
  coord_cartesian(ylim = c(0, 1))

plot raw learning rate over time for study 1

WPT1 <- WPT1 %>%
  group_by(Participant) %>%
  mutate(Trial = seq_len(n()))

#Study 1 
ggplot(WPT1, aes(Trial, acc, color = Condition)) + 
  geom_smooth(method = "loess")+
  scale_y_continuous(name = "Accuracy") 

#Study2
ggplot(WPT2, aes(Trial, acc, color = Condition)) + 
  geom_smooth(method = "loess")+
  scale_y_continuous(name = "Accuracy") 

#Study3
ggplot(WPT3, aes(Trial, acc, color = Condition)) + 
  geom_smooth(method = "loess")+
  scale_y_continuous(name = "Accuracy") 

########RL Results#######

Load RL Models

RL.ML.3 <- read.csv("fullBICStudy3.ML.csv", header=T, stringsAsFactors = FALSE, na.strings=c("","NA"))
RL.Map.3 <- read.csv("fullBICStudy3.MAP.csv", header=T, stringsAsFactors = FALSE, na.strings=c("","NA"))

RL.ML.2 <- read.csv("fullBICStudy2.ML.csv", header=T, stringsAsFactors = FALSE, na.strings=c("","NA"))
RL.Map.2 <- read.csv("fullBICStudy2.MAP.csv", header=T, stringsAsFactors = FALSE, na.strings=c("","NA"))

##Individual Differences for ET DECAY
WPT3Neg <-  read.csv("Negative.ETDecay.MAP.ParamDf.csv", header=T, stringsAsFactors = FALSE, na.strings=c("","NA"))
WPT3Neg$Cond <- "Negative"
WPT3Neg <- WPT3Neg[,-1]

WPT3Pos <-  read.csv("Positive.ETDecay.MAP.ParamDf.csv", header=T, stringsAsFactors = FALSE, na.strings=c("","NA"))
WPT3Pos$Cond <- "Positive"

WPT3FullRL <- rbind(WPT3Neg, WPT3Pos)
WPT3FullRL <- WPT3FullRL[!(WPT3FullRL$subID %in% count_NAStudy3),] #remove outliers participants

WPT2Steal <-  read.csv("Steal.ETDecay.MAP.ParamDf.csv", header=T, stringsAsFactors = FALSE, na.strings=c("","NA"))
WPT2Steal$Cond <- "steal"

WPT2Weather<-  read.csv("Weather.ETDecay.MAP.ParamDf.csv", header=T, stringsAsFactors = FALSE, na.strings=c("","NA"))
WPT2Weather$Cond <- "Weather"

WPT2StealCloud<-  read.csv("StealCl.ETDecay.MAP.ParamDf.csv", header=T, stringsAsFactors = FALSE, na.strings=c("","NA"))
WPT2StealCloud$Cond <- "StealCl"

WPT2WeatherFace<-  read.csv("WeatherFa.ETDecay.MAP.ParamDf.csv", header=T, stringsAsFactors = FALSE, na.strings=c("","NA"))
WPT2WeatherFace$Cond <- "WeatherFa"

WPT2FullRL <- rbind(WPT2Steal, WPT2Weather, WPT2StealCloud, WPT2WeatherFace)
WPT2FullRL <- WPT2FullRL[!(WPT2FullRL$subID %in% count_NAStudy2),] #remove outliers participants

Plot BIC Summed Difference From Baseline Model (i.e. model that does not take into account experimental design)

ggplot(RL.ML.2, aes(x = Cond, y = SumBIC))+
  geom_bar(stat = "identity", alpha=0.5) +facet_grid(~model)+
  geom_errorbar(aes(x=as.factor(Cond), ymin=SumBIC-SE_PL, ymax=SumBIC+SE_PL))

ggplot(RL.Map.2, aes(x = Cond, y = SumBIC))+
  geom_bar(stat = "identity", alpha=0.5) +facet_grid(~model)+
  geom_errorbar(aes(x=as.factor(Cond), ymin=SumBIC-SE_PL, ymax=SumBIC+SE_PL))

ggplot(RL.ML.3, aes(x = Cond, y = SumBIC))+
  geom_bar(stat = "identity", alpha=0.5) +facet_grid(~model)+
  geom_errorbar(aes(x=as.factor(Cond), ymin=SumBIC-SE_PL, ymax=SumBIC+SE_PL))

ggplot(RL.Map.3, aes(x = Cond, y = SumBIC))+
  geom_bar(stat = "identity", alpha=0.5) +facet_grid(~model)+
  geom_errorbar(aes(x=as.factor(Cond), ymin=SumBIC-SE_PL, ymax=SumBIC+SE_PL))

Plotting and Testing for ET Decay BIC differences

#summarize study 2 BIC 
BIC.Sum.2.ETDecay <- WPT2FullRL %>% # the names of the new data frame and the data frame to be summarised
  group_by(Cond) %>%   # the grouping variable
  dplyr::summarise(SumBIC = sum(BIC, na.rm = T),  # calculates the mean of each group
                   sd_PL = sd(BIC),
                   n_PL = n(),  # calculates the sample size per group
                   SE_PL = sd(BIC)/sqrt(n())) # calculates the standard error of each group

`summarise()` ungrouping output (override with `.groups` argument)

#summarize study 3 BIC 
BIC.Sum.3.ETDecay <- WPT3FullRL %>% # the names of the new data frame and the data frame to be summarised
  group_by(Cond) %>%   # the grouping variable
  dplyr::summarise(SumBIC = sum(BIC, na.rm = T),  # calculates the mean of each group
                   sd_PL = sd(BIC),
                   n_PL = n(),  # calculates the sample size per group
                   SE_PL = sd(BIC)/sqrt(n())) # calculates the standard error of each group

`summarise()` ungrouping output (override with `.groups` argument)

#Run ANOVA, post hoc tests and for study 2 BIC 
WPT2FullRL.AOV <- aov(WPT2FullRL$BIC~WPT2FullRL$Cond)
summary(WPT2FullRL.AOV)

                 Df  Sum Sq Mean Sq F value Pr(>F)  
WPT2FullRL$Cond   3   27858    9286   2.867 0.0365 *
Residuals       369 1195136    3239                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

TukeyHSD(WPT2FullRL.AOV)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = WPT2FullRL$BIC ~ WPT2FullRL$Cond)

$`WPT2FullRL$Cond`
                        diff       lwr       upr     p adj
StealCl-steal     -16.166597 -37.65531  5.322113 0.2125431
Weather-steal     -22.874670 -44.41835 -1.330994 0.0325028
WeatherFa-steal   -18.948126 -40.72222  2.825971 0.1130017
Weather-StealCl    -6.708074 -27.96327 14.547122 0.8476182
WeatherFa-StealCl  -2.781529 -24.27024 18.707180 0.9871347
WeatherFa-Weather   3.926544 -17.61713 25.470221 0.9655238

ggplot(BIC.Sum.2.ETDecay, aes(x=as.factor(Cond), y=SumBIC)) +
  geom_bar(stat = "identity", alpha=0.5) +
  geom_errorbar(aes(x=as.factor(Cond), ymin=SumBIC-SE_PL, ymax=SumBIC+SE_PL))+
  coord_cartesian(ylim = c(100, 25000))

#Run T.Test, post hoc tests and for study 3 BIC 
t.test(WPT3FullRL$BIC~WPT3FullRL$Cond)

    Welch Two Sample t-test

data:  WPT3FullRL$BIC by WPT3FullRL$Cond
t = 0.5537, df = 199.41, p-value = 0.5804
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -11.62906  20.70945
sample estimates:
mean in group Negative mean in group Positive 
              190.3768               185.8367 

ggplot(BIC.Sum.3.ETDecay, aes(x=as.factor(Cond), y=SumBIC)) +
  geom_bar(stat = "identity", alpha=0.5) +
  geom_errorbar(aes(x=as.factor(Cond), ymin=SumBIC-SE_PL, ymax=SumBIC+SE_PL))+
  coord_cartesian(ylim = c(100, 25000))

#Test for Differences Between Studies 2 and 3 steal conditions
Steal3 <- WPT3FullRL[which(WPT3FullRL$Cond == "Negative"),]
Steal2 <- WPT2FullRL[which(WPT2FullRL$Cond == "steal"),]
t.test(Steal3$BIC,Steal2$BIC, paired = F) # Not sig different

    Welch Two Sample t-test

data:  Steal3$BIC and Steal2$BIC
t = -0.92102, df = 191.84, p-value = 0.3582
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -24.906574   9.050359
sample estimates:
mean of x mean of y 
 190.3768  198.3050 

Plot All Correlations in stereotype congruency conditions

#Seperate Stereotype Conditions
Steal3WPT <- WPT3[which(WPT3$Condition == "neg"),]
Touchdown3 <- WPT3[which(WPT3$Condition == "pos"),]
Steal2WPT <- WPT2[which(WPT2$Condition == "steal"),]

#Isolate individual differences
colnames(Steal2WPT)

 [1] "X"               "Participant"     "LR"              "Temp"            "Anneal"          "LL"              "k"               "n"              
 [9] "BIC"             "AIC"             "Pattern"         "Wins"            "Losses"          "Super_Wins"      "Super_Losses"    "Face_Shown"     
[17] "Choice"          "RT"              "Condition"       "Pattern.Outcome" "acc"             "Choice_numeric"  "Outcome"         "Stimuli"        
[25] "Trial"           "Att"             "EX"              "HH"              "SDO"             "EMS"             "IMS"             "blk_contact"    
[33] "wht_contact"     "blk_exp"         "wht_exp"         "intergroup_anx"  "avg_acc"         "prob"            "Condition_eff"   "Cprob"          
[41] "prob_eff"       

Steal3WPT.Corr <- Steal3WPT[,c(16:30)]
Touchdown3.Corr <- Touchdown3[,c(16:30)]
Steal2WPT.Corr <- Steal2WPT[,c(3,4,5,6,27:37)]
# M <- cor(Steal3WPT.Corr, use = "pairwise.complete.obs")
# corrplot(M)
S.corr.Acc <- Steal3WPT.Corr  %>% 
  correlate()  %>% 
  focus(avg_acc)

Correlation method: 'pearson'
Missing treated using: 'pairwise.complete.obs'

S.corr.Acc %>% 
  mutate(rowname = factor(rowname, levels = rowname[order(avg_acc)])) %>%  # Order by correlation strength
  ggplot(aes(x = rowname, y = avg_acc)) +
  geom_bar(stat = "identity") +
  ylab("Correlation with Alpha") +
  xlab("Variable") + theme_grey(base_size = 8)

T.Corr.Acc <- Touchdown3.Corr  %>% 
  correlate()  %>% 
  focus(avg_acc)

Correlation method: 'pearson'
Missing treated using: 'pairwise.complete.obs'

T.Corr.Acc %>% 
  mutate(rowname = factor(rowname, levels = rowname[order(avg_acc)])) %>%  # Order by correlation strength
  ggplot(aes(x = rowname, y = avg_acc)) +
  geom_bar(stat = "identity") +
  ylab("Correlation with Alpha") +
  xlab("Variable") + theme_grey(base_size = 8)

S2.corr.Acc <- Steal2WPT.Corr  %>% 
  correlate()  %>% 
  focus(avg_acc)

Correlation method: 'pearson'
Missing treated using: 'pairwise.complete.obs'

S2.corr.Acc %>% 
  mutate(rowname = factor(rowname, levels = rowname[order(avg_acc)])) %>%  # Order by correlation strength
  ggplot(aes(x = rowname, y = avg_acc)) +
  geom_bar(stat = "identity") +
  ylab("Correlation with Alpha") +
  xlab("Variable") + theme_grey(base_size = 8)

#Scale and median split individual differences for visuals

S2 <- Steal2WPT[,c(2,5,17,18,21,22,24,25,27:38)]
S3 <- Steal3WPT[,c(2,4,5,8,9,12, 14,16:26,29,36)]
stealFull <- rbind(S2, S3)

WPT2$decay <- WPT2$Anneal
WPT2$scDecay <- scale(WPT2$decay)
WPT2$decay2 <- WPT2$scDecay^2
WPT2$SDO_sca <- scale(WPT2$SDO)
WPT2$blk_exp_sca <- scale(WPT2$blk_exp)
WPT2$blk_contact_sca <- scale(WPT2$blk_contact)
WPT2$intergroup_anx_sca <- scale(WPT2$intergroup_anx)
WPT2$IMS_sca <- scale(WPT2$IMS)
WPT2$EMS_sca <- scale(WPT2$EMS)

#scale ind diff for positive condition
WPT3$decay <- WPT3$Anneal
WPT3$SDO_sca <- scale(WPT3$SDO)
WPT3$blk_exp_sca <- scale(WPT3$blk_exp)
WPT3$blk_contact_sca <- scale(WPT3$blk_contact)
WPT3$intergroup_anx_sca <- scale(WPT3$intergroup_anx)
WPT3$IMS_sca <- scale(WPT3$IMS)
WPT3$EMS_sca <- scale(WPT3$EMS)
WPT3$scDecay <- scale(WPT3$decay)
WPT3$decay2 <- WPT3$scDecay^2

stealStudy2 <- S2
stealStudy2$decay <- stealStudy2$Anneal
stealStudy2$SDO_sca <- scale(stealStudy2$SDO)
stealStudy2$blk_exp_sca <- scale(stealStudy2$blk_exp)
stealStudy2$blk_contact_sca <- scale(stealStudy2$blk_contact)
stealStudy2$intergroup_anx_sca <- scale(stealStudy2$intergroup_anx)
stealStudy2$IMS_sca <- scale(stealStudy2$IMS)
stealStudy2$EMS_sca <- scale(stealStudy2$EMS)
stealStudy2$scDecay <- scale(stealStudy2$decay)
stealStudy2$decay2 <- stealStudy2$scDecay^2

stealStudy3 <- S3
stealStudy3$decay <- stealStudy3$Anneal
stealStudy3$SDO_sca <- scale(stealStudy3$SDO)
stealStudy3$blk_exp_sca <- scale(stealStudy3$blk_exp)
stealStudy3$blk_contact_sca <- scale(stealStudy3$blk_contact)
stealStudy3$intergroup_anx_sca <- scale(stealStudy3$intergroup_anx)
stealStudy3$IMS_sca <- scale(stealStudy3$IMS)
stealStudy3$EMS_sca <- scale(stealStudy3$EMS)
stealStudy3$scDecay <- scale(stealStudy3$decay)
stealStudy3$decay2 <- stealStudy3$scDecay^2

posStudy3 <- Touchdown3
posStudy3$decay <- posStudy3$Anneal
posStudy3$SDO_sca <- scale(posStudy3$SDO)
posStudy3$blk_exp_sca <- scale(posStudy3$blk_exp)
posStudy3$blk_contact_sca <- scale(posStudy3$blk_contact)
posStudy3$intergroup_anx_sca <- scale(posStudy3$intergroup_anx)
posStudy3$IMS_sca <- scale(posStudy3$IMS)
posStudy3$EMS_sca <- scale(posStudy3$EMS)
posStudy3$scDecay <- scale(posStudy3$decay)
posStudy3$decay2 <- posStudy3$scDecay^2

stealFull$decay <- stealFull$Anneal
stealFull$scDecay <- scale(stealFull$decay)
stealFull$decay2 <- stealFull$scDecay^2
stealFull$SDO_sca <- scale(stealFull$SDO)
stealFull$blk_exp_sca <- scale(stealFull$blk_exp)
stealFull$blk_contact_sca <- scale(stealFull$blk_contact)
stealFull$intergroup_anx_sca <- scale(stealFull$intergroup_anx)
stealFull$IMS_sca <- scale(stealFull$IMS)
stealFull$EMS_sca <- scale(stealFull$EMS)


#median split 
stealStudy2$SDODich[stealStudy2$SDO > median(stealStudy2$SDO, na.rm = T) ] <- "high"
stealStudy2$SDODich[stealStudy2$SDO < median(stealStudy2$SDO, na.rm = T) ] <- "low"
stealStudy3$SDODich[stealStudy3$SDO > median(stealStudy3$SDO, na.rm = T) ] <- "high"
stealStudy3$SDODich[stealStudy3$SDO < median(stealStudy3$SDO, na.rm = T) ] <- "low"
WPT2$SDODich[WPT2$SDO > median(stealStudy2$SDO, na.rm = T) ] <- "high"
WPT2$SDODich[WPT2$SDO < median(stealStudy2$SDO, na.rm = T) ] <- "low"
WPT3$SDODich[WPT3$SDO > median(stealStudy2$SDO, na.rm = T) ] <- "high"
WPT3$SDODich[WPT3$SDO < median(stealStudy2$SDO, na.rm = T) ] <- "low"
stealFull$SDODich[stealFull$SDO > median(stealFull$SDO, na.rm = T) ] <- "high"
stealFull$SDODich[stealFull$SDO < median(stealFull$SDO, na.rm = T) ] <- "low"
table(stealStudy2$SDODich)

high  low 
7453 7734 

table(stealStudy3$intergroup_anx)

               1              1.5 1.66666666666667 1.83333333333333                2 2.33333333333333              2.5 2.66666666666667 2.83333333333333 
            1818              606             1010              202             1212             2020             1010             1010             2424 
               3 3.16666666666667 3.33333333333333              3.5 3.66666666666667 3.83333333333333                4 
            3232              808              404             1010              404              404              808 

hist(stealStudy2$SDO)

hist(stealStudy3$intergroup_anx)

stealStudy2$IntAnxDich[stealStudy2$intergroup_anx > median(stealStudy2$intergroup_anx, na.rm = T) ] <- "high"
stealStudy2$IntAnxDich[stealStudy2$intergroup_anx < median(stealStudy2$intergroup_anx, na.rm = T) ] <- "low"
stealStudy3$IntAnxDich[stealStudy3$intergroup_anx > median(stealStudy3$intergroup_anx, na.rm = T) ] <- "high"
stealStudy3$IntAnxDich[stealStudy3$intergroup_anx < median(stealStudy3$intergroup_anx, na.rm = T) ] <- "low"
posStudy3$IntAnxDich[posStudy3$intergroup_anx > median(posStudy3$intergroup_anx, na.rm = T) ] <- "high"
posStudy3$IntAnxDich[posStudy3$intergroup_anx < median(posStudy3$intergroup_anx, na.rm = T) ] <- "low"
WPT2$IntAnxDich[WPT2$intergroup_anx > median(stealStudy3$intergroup_anx, na.rm = T) ] <- "high"
WPT2$IntAnxDich[WPT2$intergroup_anx < median(stealStudy3$intergroup_anx, na.rm = T) ] <- "low"
WPT3$IntAnxDich[WPT3$intergroup_anx > median(WPT3$intergroup_anx, na.rm = T) ] <- "high"
WPT3$IntAnxDich[WPT3$intergroup_anx < median(WPT3$intergroup_anx, na.rm = T) ] <- "low"
stealFull$IntAnxDich[stealFull$intergroup_anx > median(stealFull$intergroup_anx, na.rm = T) ] <- "high"
stealFull$IntAnxDich[stealFull$intergroup_anx < median(stealFull$intergroup_anx, na.rm = T) ] <- "low"

stealStudy2$IMSDich[stealStudy2$IMS > median(stealStudy2$IMS, na.rm = T) ] <- "high"
stealStudy2$IMSDich[stealStudy2$IMS < median(stealStudy2$IMS, na.rm = T) ] <- "low"
stealStudy3$IMSDich[stealStudy3$IMS > median(stealStudy3$IMS, na.rm = T) ] <- "high"
stealStudy3$IMSDich[stealStudy3$IMS < median(stealStudy3$IMS, na.rm = T) ] <- "low"
posStudy3$IMSDich[posStudy3$IMS > median(posStudy3$IMS, na.rm = T) ] <- "high"
posStudy3$IMSDich[posStudy3$IMS < median(posStudy3$IMS, na.rm = T) ] <- "low"
WPT2$IMSDich[WPT2$IMS > median(WPT2$IMS, na.rm = T) ] <- "high"
WPT2$IMSDich[WPT2$IMS < median(WPT2$IMS, na.rm = T) ] <- "low"
WPT3$IMSDich[WPT3$IMS > median(WPT3$IMS, na.rm = T) ] <- "high"
WPT3$IMSDich[WPT3$IMS < median(WPT3$IMS, na.rm = T) ] <- "low"
stealFull$IMSDich[stealFull$IMS > median(stealFull$IMS, na.rm = T) ] <- "high"
stealFull$IMSDich[stealFull$IMS < median(stealFull$IMS, na.rm = T) ] <- "low"

stealStudy2$EMSDich[stealStudy2$EMS > median(stealStudy2$EMS, na.rm = T) ] <- "high"
stealStudy2$EMSDich[stealStudy2$EMS < median(stealStudy2$EMS, na.rm = T) ] <- "low"
stealStudy3$EMSDich[stealStudy3$EMS >median(stealStudy3$EMS, na.rm = T) ] <- "high"
stealStudy3$EMSDich[stealStudy3$EMS < median(stealStudy3$EMS, na.rm = T) ] <- "low"
WPT2$EMSDich[WPT2$EMS > median(WPT2$EMS, na.rm = T) ] <- "high"
WPT2$EMSDich[WPT2$EMS < median(WPT2$EMS, na.rm = T) ] <- "low"
WPT3$EMSDich[WPT3$EMS > median(WPT3$EMS, na.rm = T) ] <- "high"
WPT3$EMSDich[WPT3$EMS < median(WPT3$EMS, na.rm = T) ] <- "low"
stealFull$EMSDich[stealFull$EMS > median(stealFull$EMS, na.rm = T) ] <- "high"
stealFull$EMSDich[stealFull$EMS < median(stealFull$EMS, na.rm = T) ] <- "low"

Run logistic mixed models for all three studies

#Study1
Study1Logistic<- glmer(acc~scale(Trial)*Condition_eff+ (1|Participant), data = WPT1, family = "binomial")
summary(Study1Logistic)

#Study2
Study2Logistic<- glmer(acc~scale(Trial)*Condition_eff+ (1|Participant), data = WPT2, family = "binomial")
summary(Study2Logistic)

#Study3
Pred3Full <- Pred3Full[!is.na(Pred3Full$Pattern),]
Study3Logistic<- glmer(acc~scale(Trial)*Condition_eff+ (1|Participant), data = WPT3, family = "binomial")
summary(Study3Logistic)

#individual differences

StealFullEMS <- glmer(acc~scale(Trial)*EMS_sca + (1|Participant), data = stealFull, family = "binomial")
summary(StealFullEMS)

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: binomial  ( logit )
Formula: acc ~ scale(Trial) * EMS_sca + (1 | Participant)
   Data: stealFull

     AIC      BIC   logLik deviance df.resid 
 33708.3  33750.6 -16849.1  33698.3    34791 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-5.3411  0.2232  0.3867  0.5481  1.4412 

Random effects:
 Groups      Name        Variance Std.Dev.
 Participant (Intercept) 0.4995   0.7067  
Number of obs: 34796, groups:  Participant, 176

Fixed effects:
                     Estimate Std. Error z value Pr(>|z|)    
(Intercept)           1.42656    0.05536  25.769  < 2e-16 ***
scale(Trial)          0.48908    0.01425  34.316  < 2e-16 ***
EMS_sca              -0.05797    0.05159  -1.123    0.261    
scale(Trial):EMS_sca -0.06522    0.01408  -4.631 3.64e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) scl(T) EMS_sc
scale(Tril)  0.066              
EMS_sca     -0.007 -0.013       
scl(T):EMS_ -0.012 -0.042  0.071

StealFullIMS <- glmer(acc~scale(Trial)*IMS_sca + (1|Participant), data = stealFull, family = "binomial")
summary(StealFullIMS)

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: binomial  ( logit )
Formula: acc ~ scale(Trial) * IMS_sca + (1 | Participant)
   Data: stealFull

     AIC      BIC   logLik deviance df.resid 
 33700.0  33742.3 -16845.0  33690.0    34793 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-5.6722  0.2210  0.3858  0.5510  1.5011 

Random effects:
 Groups      Name        Variance Std.Dev.
 Participant (Intercept) 0.5017   0.7083  
Number of obs: 34798, groups:  Participant, 176

Fixed effects:
                     Estimate Std. Error z value Pr(>|z|)    
(Intercept)           1.42866    0.05549  25.746  < 2e-16 ***
scale(Trial)          0.48926    0.01428  34.260  < 2e-16 ***
IMS_sca              -0.12278    0.05467  -2.246 0.024718 *  
scale(Trial):IMS_sca -0.05301    0.01407  -3.766 0.000166 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) scl(T) IMS_sc
scale(Tril)  0.067              
IMS_sca     -0.013 -0.014       
scl(T):IMS_ -0.014 -0.079  0.069

StealFullIntAnx <- glmer(acc~scale(Trial)*intergroup_anx_sca + (1|Participant), data = stealFull, family = "binomial")
summary(StealFullIntAnx)

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: binomial  ( logit )
Formula: acc ~ scale(Trial) * intergroup_anx_sca + (1 | Participant)
   Data: stealFull

     AIC      BIC   logLik deviance df.resid 
 33487.0  33529.3 -16738.5  33477.0    34806 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-5.7075  0.2266  0.3845  0.5439  1.3748 

Random effects:
 Groups      Name        Variance Std.Dev.
 Participant (Intercept) 0.506    0.7113  
Number of obs: 34811, groups:  Participant, 176

Fixed effects:
                                Estimate Std. Error z value Pr(>|z|)    
(Intercept)                      1.44585    0.05572  25.950   <2e-16 ***
scale(Trial)                     0.48234    0.01430  33.735   <2e-16 ***
intergroup_anx_sca              -0.06223    0.05003  -1.244    0.214    
scale(Trial):intergroup_anx_sca -0.00179    0.01444  -0.124    0.901    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) scl(T) intr__
scale(Tril)  0.066              
intrgrp_nx_ -0.007 -0.005       
scl(Trl):__ -0.006 -0.041  0.077

StealFullblkCon <- glmer(acc~scale(Trial)*blk_contact_sca + (1|Participant), data = stealFull, family = "binomial")
summary(StealFullblkCon)

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: binomial  ( logit )
Formula: acc ~ scale(Trial) * blk_contact_sca + (1 | Participant)
   Data: stealFull

     AIC      BIC   logLik deviance df.resid 
 34275.4  34317.8 -17132.7  34265.4    35387 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-5.7859  0.2238  0.3866  0.5485  1.3902 

Random effects:
 Groups      Name        Variance Std.Dev.
 Participant (Intercept) 0.514    0.7169  
Number of obs: 35392, groups:  Participant, 179

Fixed effects:
                             Estimate Std. Error z value Pr(>|z|)    
(Intercept)                   1.42795    0.05563  25.668   <2e-16 ***
scale(Trial)                  0.48111    0.01410  34.120   <2e-16 ***
blk_contact_sca              -0.01249    0.05541  -0.225    0.822    
scale(Trial):blk_contact_sca -0.01617    0.01398  -1.157    0.247    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) scl(T) blk_c_
scale(Tril)  0.064              
blk_cntct_s -0.001 -0.002       
scl(Trl):__ -0.002 -0.016  0.059

StealFullblkExt <- glmer(acc~scale(Trial)*blk_exp_sca + (1|Participant), data = stealFull, family = "binomial")
summary(StealFullblkExt)

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: binomial  ( logit )
Formula: acc ~ scale(Trial) * blk_exp_sca + (1 | Participant)
   Data: stealFull

     AIC      BIC   logLik deviance df.resid 
 31791.2  31833.2 -15890.6  31781.2    32844 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-5.7896  0.2246  0.3874  0.5470  1.4190 

Random effects:
 Groups      Name        Variance Std.Dev.
 Participant (Intercept) 0.4946   0.7033  
Number of obs: 32849, groups:  Participant, 166

Fixed effects:
                         Estimate Std. Error z value Pr(>|z|)    
(Intercept)               1.42843    0.05674  25.176   <2e-16 ***
scale(Trial)              0.48953    0.01467  33.378   <2e-16 ***
blk_exp_sca              -0.06330    0.05640  -1.122   0.2617    
scale(Trial):blk_exp_sca -0.03548    0.01439  -2.466   0.0136 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) scl(T) blk_x_
scale(Tril)  0.066              
blk_exp_sca -0.001 -0.008       
scl(Trl):__ -0.008 -0.041  0.063

RL differences by study & condition

Study2RLRT <- lmer(RT~scDecay*Condition_eff +(1|Participant), data = WPT2)
Study2RLRT <- lmer(RT~scDecay*Condition_eff +(1|Participant), data = WPT2)
summary(Study2RLRT)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: RT ~ scDecay * Condition_eff + (1 | Participant)
   Data: WPT2

REML criterion at convergence: 323716.5

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-1.380 -0.293 -0.120  0.151 78.768 

Random effects:
 Groups      Name        Variance Std.Dev.
 Participant (Intercept) 0.1327   0.3643  
 Residual                4.7471   2.1788  
Number of obs: 73482, groups:  Participant, 373

Fixed effects:
                                    Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)                          1.93603    0.02085 363.87350  92.849   <2e-16 ***
scDecay                              0.02192    0.02102 364.05420   1.043   0.2977    
Condition_effsteal                   0.04931    0.03642 363.73749   1.354   0.1766    
Condition_effsteal_clouds            0.07717    0.03651 363.92471   2.113   0.0352 *  
Condition_effweather_faces          -0.06511    0.03571 363.91190  -1.823   0.0691 .  
scDecay:Condition_effsteal          -0.06178    0.03565 364.38103  -1.733   0.0839 .  
scDecay:Condition_effsteal_clouds    0.03091    0.03851 364.31145   0.803   0.4227    
scDecay:Condition_effweather_faces  -0.02623    0.03721 363.74537  -0.705   0.4812    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
                 (Intr) scDecy Cndtn_ Cndtn_ffs_ Cndtn_ffw_ scD:C_ scDcy:Cndtn_ffs_
scDecay           0.017                                                            
Cndtn_ffstl       0.015 -0.113                                                     
Cndtn_ffst_       0.019  0.118 -0.346                                              
Cndtn_ffwt_      -0.020  0.088 -0.332 -0.333                                       
scDcy:Cndt_      -0.116 -0.036 -0.115 -0.009      0.010                            
scDcy:Cndtn_ffs_  0.113  0.098 -0.008  0.145     -0.120     -0.359                 
scDcy:Cndtn_ffw_  0.085  0.038  0.009 -0.121      0.117     -0.335 -0.384          

Plot individual differences

stealStudy2$EMSDich

   [1] "low"  "high" "low"  "low"  "low"  "low"  "low"  "high" "high" NA     NA     "low"  "low"  "low"  "high" "low"  "low"  "high" "high" "low"  "low" 
  [22] "high" "low"  "high" "low"  "low"  "high" "high" "low"  "low"  "high" "low"  "low"  "low"  NA     NA     "high" "high" "low"  "high" "low"  NA    
  [43] "high" "high" "low"  "high" "high" "high" "high" "high" "low"  "low"  "high" "high" "high" "low"  "high" "low"  "low"  "high" "low"  "low"  "low" 
  [64] "low"  NA     "low"  "high" NA     "high" NA     "high" "high" "low"  "low"  "high" "high" "low"  "high" "high" "low"  "low"  "low"  "high" "low" 
  [85] "high" "high" "low"  "high" "low"  "low"  "low"  "low"  "low"  "high" "high" NA     NA     "low"  "low"  "low"  "high" "low"  "low"  "high" "high"
 [106] "low"  "low"  "high" "low"  "high" "low"  "low"  "high" "high" "low"  "low"  "high" "low"  "low"  "low"  NA     NA     "high" "high" "low"  "high"
 [127] "low"  NA     "high" "high" "low"  "high" "high" "high" "high" "high" "low"  "low"  "high" "high" "high" "low"  "high" "low"  "low"  "high" "low" 
 [148] "low"  "low"  "low"  NA     "low"  "high" NA     "high" NA     "high" "high" "low"  "low"  "high" "high" "low"  "high" "high" "low"  "low"  "low" 
 [169] "high" "low"  "high" "high" "low"  "high" "low"  "low"  "low"  "low"  "low"  "high" "high" NA     NA     "low"  "low"  "low"  "high" "low"  "low" 
 [190] "high" "high" "low"  "low"  "high" "low"  "high" "low"  "low"  "high" "high" "low"  "low"  "high" "low"  "low"  "low"  NA     NA     "high" "high"
 [211] "low"  "high" "low"  NA     "high" "high" "low"  "high" "high" "high" "high" "high" "low"  "low"  "high" "high" "high" "low"  "high" "low"  "low" 
 [232] "high" "low"  "low"  "low"  "low"  NA     "low"  "high" NA     "high" NA     "high" "high" "low"  "low"  "high" "high" "low"  "high" "high" "low" 
 [253] "low"  "low"  "high" "low"  "high" "high" "low"  "high" "low"  "low"  "low"  "low"  "low"  "high" "high" NA     NA     "low"  "low"  "low"  "high"
 [274] "low"  "low"  "high" "high" "low"  "low"  "high" "low"  "high" "low"  "low"  "high" "high" "low"  "low"  "high" "low"  "low"  "low"  NA     NA    
 [295] "high" "high" "low"  "high" "low"  NA     "high" "high" "low"  "high" "high" "high" "high" "high" "low"  "low"  "high" "high" "high" "low"  "high"
 [316] "low"  "low"  "high" "low"  "low"  "low"  "low"  NA     "low"  "high" NA     "high" NA     "high" "high" "low"  "low"  "high" "high" "low"  "high"
 [337] "high" "low"  "low"  "low"  "high" "low"  "high" "high" "low"  "high" "low"  "low"  "low"  "low"  "low"  "high" "high" NA     NA     "low"  "low" 
 [358] "low"  "high" "low"  "low"  "high" "high" "low"  "low"  "high" "low"  "high" "low"  "low"  "high" "high" "low"  "low"  "high" "low"  "low"  "low" 
 [379] NA     NA     "high" "high" "low"  "high" "low"  NA     "high" "high" "low"  "high" "high" "high" "high" "high" "low"  "low"  "high" "high" "high"
 [400] "low"  "high" "low"  "low"  "high" "low"  "low"  "low"  "low"  NA     "low"  "high" NA     "high" NA     "high" "high" "low"  "low"  "high" "high"
 [421] "low"  "high" "high" "low"  "low"  "low"  "high" "low"  "high" "high" "low"  "high" "low"  "low"  "low"  "low"  "low"  "high" "high" NA     NA    
 [442] "low"  "low"  "low"  "high" "low"  "low"  "high" "high" "low"  "low"  "high" "low"  "high" "low"  "low"  "high" "high" "low"  "low"  "high" "low" 
 [463] "low"  "low"  NA     NA     "high" "high" "low"  "high" "low"  NA     "high" "high" "low"  "high" "high" "high" "high" "high" "low"  "low"  "high"
 [484] "high" "high" "low"  "high" "low"  "low"  "high" "low"  "low"  "low"  "low"  NA     "low"  "high" NA     "high" NA     "high" "high" "low"  "low" 
 [505] "high" "high" "low"  "high" "high" "low"  "low"  "low"  "high" "low"  "high" "high" "low"  "high" "low"  "low"  "low"  "low"  "low"  "high" "high"
 [526] NA     NA     "low"  "low"  "low"  "high" "low"  "low"  "high" "high" "low"  "low"  "high" "low"  "high" "low"  "low"  "high" "high" "low"  "low" 
 [547] "high" "low"  "low"  "low"  NA     NA     "high" "high" "low"  "high" "low"  NA     "high" "high" "low"  "high" "high" "high" "high" "high" "low" 
 [568] "low"  "high" "high" "high" "low"  "high" "low"  "low"  "high" "low"  "low"  "low"  "low"  NA     "low"  "high" NA     "high" NA     "high" "high"
 [589] "low"  "low"  "high" "high" "low"  "high" "high" "low"  "low"  "low"  "high" "low"  "high" "high" "low"  "high" "low"  "low"  "low"  "low"  "low" 
 [610] "high" "high" NA     NA     "low"  "low"  "low"  "high" "low"  "low"  "high" "high" "low"  "low"  "high" "low"  "high" "low"  "low"  "high" "high"
 [631] "low"  "low"  "high" "low"  "low"  "low"  NA     NA     "high" "high" "low"  "high" "low"  NA     "high" "high" "low"  "high" "high" "high" "high"
 [652] "high" "low"  "low"  "high" "high" "high" "low"  "high" "low"  "low"  "high" "low"  "low"  "low"  "low"  NA     "low"  "high" NA     "high" NA    
 [673] "high" "high" "low"  "low"  "high" "high" "low"  "high" "high" "low"  "low"  "low"  "high" "low"  "high" "high" "low"  "high" "low"  "low"  "low" 
 [694] "low"  "low"  "high" "high" NA     NA     "low"  "low"  "low"  "high" "low"  "low"  "high" "high" "low"  "low"  "high" "low"  "high" "low"  "low" 
 [715] "high" "high" "low"  "low"  "high" "low"  "low"  "low"  NA     NA     "high" "high" "low"  "high" "low"  NA     "high" "high" "low"  "high" "high"
 [736] "high" "high" "high" "low"  "low"  "high" "high" "high" "low"  "high" "low"  "low"  "high" "low"  "low"  "low"  "low"  NA     "low"  "high" NA    
 [757] "high" NA     "high" "high" "low"  "low"  "high" "high" "low"  "high" "high" "low"  "low"  "low"  "high" "low"  "high" "high" "low"  "high" "low" 
 [778] "low"  "low"  "low"  "low"  "high" "high" NA     NA     "low"  "low"  "low"  "high" "low"  "low"  "high" "high" "low"  "low"  "high" "low"  "high"
 [799] "low"  "low"  "high" "high" "low"  "low"  "high" "low"  "low"  "low"  NA     NA     "high" "high" "low"  "high" "low"  NA     "high" "high" "low" 
 [820] "high" "high" "high" "high" "high" "low"  "low"  "high" "high" "high" "low"  "high" "low"  "low"  "high" "low"  "low"  "low"  "low"  NA     "low" 
 [841] "high" NA     "high" NA     "high" "high" "low"  "low"  "high" "high" "low"  "high" "high" "low"  "low"  "low"  "high" "low"  "high" "high" "low" 
 [862] "high" "low"  "low"  "low"  "low"  "low"  "high" "high" NA     NA     "low"  "low"  "low"  "high" "low"  "low"  "high" "high" "low"  "low"  "high"
 [883] "low"  "high" "low"  "low"  "high" "high" "low"  "low"  "high" "low"  "low"  "low"  NA     NA     "high" "high" "low"  "high" "low"  NA     "high"
 [904] "high" "low"  "high" "high" "high" "high" "high" "low"  "low"  "high" "high" "high" "low"  "high" "low"  "low"  "high" "low"  "low"  "low"  "low" 
 [925] NA     "low"  "high" NA     "high" NA     "high" "high" "low"  "low"  "high" "high" "low"  "high" "high" "low"  "low"  "low"  "high" "low"  "high"
 [946] "high" "low"  "high" "low"  "low"  "low"  "low"  "low"  "high" "high" NA     NA     "low"  "low"  "low"  "high" "low"  "low"  "high" "high" "low" 
 [967] "low"  "high" "low"  "high" "low"  "low"  "high" "high" "low"  "low"  "high" "low"  "low"  "low"  NA     NA     "high" "high" "low"  "high" "low" 
 [988] NA     "high" "high" "low"  "high" "high" "high" "high" "high" "low"  "low"  "high" "high"
 [ reached getOption("max.print") -- omitted 15232 entries ]