Accuracy Analyses for the Discrimination Session
Data Files
rm(list=ls())
#libraries
library(dplyr)
library(lattice)
library(languageR)
library(mgcv)
library(itsadug)
library(parallel)
# All trials in data, apart from practice trials
# Load Data Files
# Category Learning
data1<-read.table("TrialRep_StimPresentationA2_CAT.txt", header=T)
# Paired-Associate Learning
data2<-read.table("TrialRep_StimPresentationA2_PA.txt", header=T)
# Combine data files
data<-rbind(data1, data2)
str(data)
## 'data.frame': 4608 obs. of 17 variables:
## $ Sbj_name : chr "aggfyt95" "aggfyt95" "aggfyt95" "aggfyt95" ...
## $ Group : chr "c" "c" "c" "c" ...
## $ subcondition : chr "s4" "s4" "s4" "s4" ...
## $ IP_LABEL : chr "StimPresentation" "StimPresentation" "StimPresentation" "StimPresentation" ...
## $ Trial_Number : int 25 26 27 28 29 30 31 32 33 34 ...
## $ condition : chr "mixed" "mixed" "ideog" "label" ...
## $ same_diff : chr "diff" "diff" "same" "same" ...
## $ cond_left : chr "ideog" "ideog" "ideog" "label" ...
## $ cond_right : chr "label" "label" "ideog" "label" ...
## $ Response : chr "Correct" "Correct" "Correct" "Correct" ...
## $ BLINK_COUNT : int 1 1 1 1 1 1 1 1 1 2 ...
## $ FIXATION_COUNT : int 7 6 9 9 11 9 10 9 6 9 ...
## $ FIXATION_DURATION_MIN : chr "87" "152" "110" "115" ...
## $ FIXATION_DURATION_MAX : chr "756" "1264" "506" "459" ...
## $ RUN_COUNT : int 4 3 6 6 5 4 6 5 3 4 ...
## $ VISITED_INTEREST_AREA_COUNT: int 2 2 2 2 2 2 2 2 2 2 ...
## $ RT : int 1316 1354 1347 1495 1074 936 1436 1428 1343 1355 ...
# Convert all chr columns to factor:
data <- mutate_if(data, is.character, as.factor)
# Rename Group levels
levels(data$Group)<-c("CAT", "PA")
Preprocessing
# In Experiment Builder, no response is denoted by a negative RT. Make necessary conversions
data$RT<-ifelse(data$RT>0, data$RT, NA)
sum(is.na(data$RT))
## [1] 16
# Adjust 'Response' Column
data$Response<-as.factor(ifelse(is.na(data$RT),NA,data$Response))
sum(is.na(data$Response))
## [1] 16
levels(data$Response)
## [1] "1" "2"
levels(data$Response)<-c('Correct', 'Wrong')
sum(is.na(data$Response))
## [1] 16
# Create 'acc' Variable, rename Variables
data$acc<-ifelse(is.na(data$Response), NA, ifelse(data$Response=='Correct', 1,0))
colnames(data)[1] <- "sbj"
colnames(data)[6] <- "cnd"
#################################
# keep correct responses only #
#################################
sum(ifelse(data$acc %in% 0, 1,0))
## [1] 83
sum(ifelse(data$acc %in% 0 & data$Group=='CAT', 1,0))
## [1] 43
sum(ifelse(data$acc %in% 0 & data$Group=='PA', 1,0))
## [1] 40
83/4608 #percent of trials
## [1] 0.01801215
sum(ifelse(is.na(data$Response) | data$Response %in% 'Wrong', 1,0))
## [1] 99
#total number of data points excluded.
data_cr<-droplevels(subset(data, Response == 'Correct'))
# LogRT
data_cr$logRT<-log(data_cr$RT)
# Create trial variable
data_cr$trial<-data_cr$Trial_Number-24 #first 24 trials are practice trials, removed from the report
range(data_cr$trial, na.rm=T)
## [1] 1 96
sum(is.na(data_cr$trial))
## [1] 0
str(data_cr)
## 'data.frame': 4509 obs. of 20 variables:
## $ sbj : Factor w/ 48 levels "aggfyt95","agkart97",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Group : Factor w/ 2 levels "CAT","PA": 1 1 1 1 1 1 1 1 1 1 ...
## $ subcondition : Factor w/ 8 levels "s1","s2","s3",..: 4 4 4 4 4 4 4 4 4 4 ...
## $ IP_LABEL : Factor w/ 1 level "StimPresentation": 1 1 1 1 1 1 1 1 1 1 ...
## $ Trial_Number : int 25 26 27 28 29 30 32 33 34 35 ...
## $ cnd : Factor w/ 3 levels "ideog","label",..: 3 3 1 2 1 3 1 2 1 3 ...
## $ same_diff : Factor w/ 2 levels "diff","same": 1 1 2 2 2 1 1 1 2 1 ...
## $ cond_left : Factor w/ 2 levels "ideog","label": 1 1 1 2 1 1 1 2 1 2 ...
## $ cond_right : Factor w/ 2 levels "ideog","label": 2 2 1 2 1 2 1 2 1 1 ...
## $ Response : Factor w/ 1 level "Correct": 1 1 1 1 1 1 1 1 1 1 ...
## $ BLINK_COUNT : int 1 1 1 1 1 1 1 1 2 1 ...
## $ FIXATION_COUNT : int 7 6 9 9 11 9 9 6 9 4 ...
## $ FIXATION_DURATION_MIN : Factor w/ 531 levels "0","0,00","1",..: 503 126 31 41 99 95 164 65 519 194 ...
## $ FIXATION_DURATION_MAX : Factor w/ 2028 levels "1000","1000,00",..: 1633 308 1156 1063 841 1748 2026 1953 800 544 ...
## $ RUN_COUNT : int 4 3 6 6 5 4 5 3 4 3 ...
## $ VISITED_INTEREST_AREA_COUNT: int 2 2 2 2 2 2 2 2 2 2 ...
## $ RT : int 1316 1354 1347 1495 1074 936 1428 1343 1355 1434 ...
## $ acc : num 1 1 1 1 1 1 1 1 1 1 ...
## $ logRT : num 7.18 7.21 7.21 7.31 6.98 ...
## $ trial : num 1 2 3 4 5 6 8 9 10 11 ...
xylowess.fnc(RT~trial|sbj, data=data_cr)
#It seems that RTs are dependent on trial number
#Previous research suggests a difference between same and different trials
#We thus analyse RT data using same/different and also Group (Category vs. Paired-associate) as predictors
Graph
#############################################
#Calculate average RT by Group and Same_Diff
#############################################
data_rt<-aggregate(RT~Group+same_diff+trial, data_cr, mean, na.rm=T)
color=c("red", "#D0A0A0","blue","cornflowerblue")
lt=c(1,1,2,2)
rng<-range(data_rt$RT)
#rng[2] is 1492.818 thus, ylim is set to 1500 ms
#CAT
plot(RT~trial, data_rt, subset=(same_diff=="same"&Group=="CAT"), type="l", main="Discrimination, Mean Response Time\nby Group and Same/Different", ylab="RT (msec)", xlab="Trial", col=color[1], frame.plot=F, ylim=c(0,1500), las=2, cex.axis=0.8, lwd=2, lty=lt[1], xaxt="n")
lines(RT~jitter(trial), data_rt, subset=(same_diff=="diff"&Group=="CAT"), col=color[2], lwd=2, lty=lt[2])
#PA
lines(RT~jitter(trial), data_rt, subset=(same_diff=="same"&Group=="PA"), col=color[3], lwd=2, lty=lt[3])
lines(RT~jitter(trial), data_rt, subset=(same_diff=="diff"&Group=="PA"), col=color[4], lwd=2, lty=lt[4])
axis(side=1, at = c(1,24,48,72,96), labels = c(1,24,48,64,96), cex.axis = 0.8)
legend(x=15, y=500, legend=c("Category Learning, Same", "Category Learning, Different","Paired-Associate, Same","Paired-Associate, Different"), lty=lt, bty="n", col=color, seg.len=4.0, cex=1)
#########################
RTs by Learning Group
data_group0<-aggregate(RT~Group+sbj, data_cr, mean, na.rm=T)
data_groupmean<-aggregate(RT~Group, data_group0, mean, na.rm=T)
data_groupsd<-aggregate(RT~Group, data_group0, se, na.rm=T)
t.test(RT~Group, data=data_group0)
##
## Welch Two Sample t-test
##
## data: RT by Group
## t = 2.2625, df = 41.361, p-value = 0.02898
## alternative hypothesis: true difference in means between group CAT and group PA is not equal to 0
## 95 percent confidence interval:
## 12.81797 225.36212
## sample estimates:
## mean in group CAT mean in group PA
## 1081.7109 962.6209
#there's a difference between groups
RTs by Learning Group and Same/Different
#Category-learning group
datagroupsamecat<-aggregate(RT~same_diff+sbj, data_cr[data_cr$Group=="CAT",], mean, na.rm=T)
t.test(RT~same_diff, data=datagroupsamecat)
##
## Welch Two Sample t-test
##
## data: RT by same_diff
## t = -0.86757, df = 43.865, p-value = 0.3904
## alternative hypothesis: true difference in means between group diff and group same is not equal to 0
## 95 percent confidence interval:
## -188.8216 75.1835
## sample estimates:
## mean in group diff mean in group same
## 1053.932 1110.751
#no difference
#Paired-Associate group
datagroupsamepa<-aggregate(RT~same_diff+sbj, data_cr[data_cr$Group=="PA",], mean, na.rm=T)
t.test(RT~same_diff, data=datagroupsamepa)
##
## Welch Two Sample t-test
##
## data: RT by same_diff
## t = -0.20498, df = 44.865, p-value = 0.8385
## alternative hypothesis: true difference in means between group diff and group same is not equal to 0
## 95 percent confidence interval:
## -99.60594 81.20572
## sample estimates:
## mean in group diff mean in group same
## 958.3291 967.5292
#no difference
Main Analysis
library(fitdistrplus)
descdist(data_cr$RT, boot = 1000)
## summary statistics
## ------
## min: 463 max: 2926
## median: 939
## mean: 1020.595
## estimated sd: 339.9001
## estimated skewness: 1.756181
## estimated kurtosis: 7.56242
fg <- fitdist(data_cr$RT, "gamma")
fln <- fitdist(data_cr$RT, "lnorm")
#dev.off()
par(mfrow = c(2, 2))
plot.legend <- c( "lognormal", "gamma")
denscomp(list( fln, fg), legendtext = plot.legend)
qqcomp(list( fln, fg), legendtext = plot.legend)
cdfcomp(list( fln, fg), legendtext = plot.legend)
ppcomp(list( fln, fg), legendtext = plot.legend)
par(mfrow = c(1, 1))
#dev.off()
library(lme4)
m0 = glmer(RT~ (1|sbj),
family ="Gamma"(link=identity),data =data_cr )
m1 = glmer(RT~ Group+(1|sbj),
family ="Gamma"(link=identity),data =data_cr )
anova(m0,m1)
## Data: data_cr
## Models:
## m0: RT ~ (1 | sbj)
## m1: RT ~ Group + (1 | sbj)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m0 3 62340 62359 -31167 62334
## m1 4 62338 62364 -31165 62330 3.8272 1 0.05043 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#m0 preferred
m2 = glmer(RT~ Group+same_diff+(1|sbj),
family ="Gamma"(link=identity),data =data_cr )
anova(m0,m2)
## Data: data_cr
## Models:
## m0: RT ~ (1 | sbj)
## m2: RT ~ Group + same_diff + (1 | sbj)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m0 3 62340 62359 -31167 62334
## m2 5 62327 62359 -31158 62317 17.244 2 0.0001801 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#m2 preferred
m3 = glmer(RT~ trial*Group+same_diff+(1|sbj),
family ="Gamma"(link=identity),data =data_cr )
anova(m2,m3)
## Data: data_cr
## Models:
## m2: RT ~ Group + same_diff + (1 | sbj)
## m3: RT ~ trial * Group + same_diff + (1 | sbj)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m2 5 62327 62359 -31158 62317
## m3 7 62215 62260 -31101 62201 115.46 2 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#m3 preferred
m4= glmer(RT~ trial*Group+trial*same_diff+(1|sbj),
family ="Gamma"(link=identity),data =data_cr )
anova(m3,m4)
## Data: data_cr
## Models:
## m3: RT ~ trial * Group + same_diff + (1 | sbj)
## m4: RT ~ trial * Group + trial * same_diff + (1 | sbj)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m3 7 62215 62260 -31101 62201
## m4 8 62193 62245 -31089 62177 23.885 1 1.023e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#m4 preferred
m5<-glmer(RT~ trial*Group*same_diff+(1|sbj),
family ="Gamma"(link=identity),data =data_cr )
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0262156 (tol = 0.002, component 1)
#failure to converge
fm<-m4
#final model
summary(fm)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Gamma ( identity )
## Formula: RT ~ trial * Group + trial * same_diff + (1 | sbj)
## Data: data_cr
##
## AIC BIC logLik deviance df.resid
## 62193.5 62244.8 -31088.7 62177.5 4501
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8938 -0.6284 -0.2041 0.3867 10.0871
##
## Random effects:
## Groups Name Variance Std.Dev.
## sbj (Intercept) 8.107e+03 90.0409
## Residual 7.106e-02 0.2666
## Number of obs: 4509, groups: sbj, 48
##
## Fixed effects:
## Estimate Std. Error t value Pr(>|z|)
## (Intercept) 1151.3789 7.2965 157.800 < 2e-16 ***
## trial -1.0317 0.1944 -5.307 1.11e-07 ***
## GroupPA -159.3144 5.7487 -27.713 < 2e-16 ***
## same_diffsame 84.4108 5.2904 15.955 < 2e-16 ***
## trial:GroupPA 0.5839 0.2435 2.398 0.0165 *
## trial:same_diffsame -1.2082 0.1438 -8.401 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trial GropPA sm_dff tr:GPA
## trial 0.003
## GroupPA 0.096 0.011
## same_diffsm 0.167 0.187 0.029
## trial:GrpPA 0.002 -0.687 -0.026 -0.002
## trl:sm_dffs -0.088 -0.377 -0.015 -0.554 0.008
library(sjPlot)
#install.packages("Rtools")
levels(data_cr$Group)
## [1] "CAT" "PA"
levels(data_cr$same_diff)
## [1] "diff" "same"
#CAT vs. PA for "different" trials
plot_model(fm, type=c("int"))
## [[1]]
##
## [[2]]
#same vs. different trials for the Category-learning group
# relevel the Group factor to get same vs. different for the Paired-Associate group
data_cr$Group<-factor(data_cr$Group, levels=c("PA", "CAT"))
levels(data_cr$Group)
## [1] "PA" "CAT"
fm_= glmer(RT~ trial*Group+trial*same_diff+(1|sbj),
family ="Gamma"(link=identity),data =data_cr, control=glmerControl(optCtrl=list(maxfun=3e5)) )
#failure to converge
# relevel the same_diff factor (and Group as it was) to get CAT vs, PA for "same" trials
data_cr$Group<-factor(data_cr$Group, levels=c("CAT", "PA"))
data_cr$same_diff<-factor(data_cr$same_diff, levels=c("same", "diff"))
levels(data_cr$same_diff)
## [1] "same" "diff"
fm__<-glmer(RT~ trial*Group+trial*same_diff+(1|sbj),
family ="Gamma"(link=identity),data =data_cr, control=glmerControl(optCtrl=list(maxfun=3e5)) )
# CAT vs. PA for "same" trials
plot_model(fm__, type=c("int"))
## [[1]]
##
## [[2]]
#same vs. different trials for the Category-learning group (same as previous)
# Visualisation of model estimates suggest that RTs depend on trial.
# More specifically, the speed of responding increased as trials progressed.
Session Information
sessionInfo()
## R version 4.2.1 (2022-06-23 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19044)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=Greek_Greece.utf8 LC_CTYPE=Greek_Greece.utf8
## [3] LC_MONETARY=Greek_Greece.utf8 LC_NUMERIC=C
## [5] LC_TIME=Greek_Greece.utf8
##
## attached base packages:
## [1] parallel stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] sjPlot_2.8.10 lme4_1.1-29 Matrix_1.4-1 fitdistrplus_1.1-8
## [5] survival_3.3-1 MASS_7.3-57 itsadug_2.4.1 plotfunctions_1.4
## [9] mgcv_1.8-40 nlme_3.1-157 languageR_1.5.0 lattice_0.20-45
## [13] dplyr_1.0.9
##
## loaded via a namespace (and not attached):
## [1] Rcpp_1.0.8.3 mvtnorm_1.1-3 tidyr_1.2.0 digest_0.6.29
## [5] utf8_1.2.2 R6_2.5.1 backports_1.4.1 evaluate_0.15
## [9] coda_0.19-4 ggplot2_3.3.6 highr_0.9 pillar_1.7.0
## [13] rlang_1.0.2 rstudioapi_0.13 minqa_1.2.4 performance_0.9.1
## [17] nloptr_2.0.3 jquerylib_0.1.4 effectsize_0.7.0 rmarkdown_2.14
## [21] labeling_0.4.2 ggeffects_1.1.2 splines_4.2.1 stringr_1.4.0
## [25] munsell_0.5.0 broom_0.8.0 modelr_0.1.8 compiler_4.2.1
## [29] xfun_0.31 pkgconfig_2.0.3 parameters_0.18.1 htmltools_0.5.2
## [33] insight_0.17.1 tidyselect_1.1.2 tibble_3.1.7 fansi_1.0.3
## [37] crayon_1.5.1 sjmisc_2.8.9 grid_4.2.1 jsonlite_1.8.0
## [41] xtable_1.8-4 gtable_0.3.0 lifecycle_1.0.1 magrittr_2.0.3
## [45] bayestestR_0.12.1 scales_1.2.0 datawizard_0.4.1 estimability_1.3
## [49] cli_3.3.0 stringi_1.7.6 farver_2.1.0 bslib_0.3.1
## [53] ellipsis_0.3.2 generics_0.1.2 vctrs_0.4.1 sjlabelled_1.2.0
## [57] boot_1.3-28 RColorBrewer_1.1-3 tools_4.2.1 glue_1.6.2
## [61] purrr_0.3.4 sjstats_0.18.1 emmeans_1.7.5 fastmap_1.1.0
## [65] yaml_2.3.5 colorspace_2.0-3 knitr_1.39 sass_0.4.1