### Workspace ------------------------------------------------------------------- # To reset the workspace rm(list=ls()) ### Loading packages ------------------------------------------------------------ library(plyr) library(dplyr) library(gdata) library(reshape) library(car) library(tidyr) library(magrittr) library(data.table) library(ggplot2) # for the graphics library(readr) library(stringr) library(psych) library(GPArotation) # Loading a command for the confidence intervals ci.fun <- function (x) { moe<-qt(0.975, length(x)-1) * sd(x,na.rm=TRUE) / sqrt(length(x)) LI <- mean(x,na.rm=TRUE) - moe HI <- mean(x,na.rm=TRUE) + moe Mean <- mean(x,na.rm=TRUE) CI<-cbind(LI,Mean,HI) colnames(CI)<-c("2.5%","Mean","97.5%") row.names(CI)<-deparse(substitute(x)) CI } # Ensure that relative paths start from the same directory as this script rstudioapi::getActiveDocumentContext()$path %>% dirname %>% setwd ### Data files import and dataset formating ----------------------------------------------------------- # Loading the Rdata file (created from the readPsytoolkit.R file) load("dataVAAST.RData") DF <- VAAST # Recoding the variables modalities with a clearer signification DF$Movement <- recode (DF$Movement, "'1'='avoid'; '2'='appr'") DF$Category <- recode (DF$Category, "'1'='catA'; '2'='catB'") DF$BlockOrder <- recode (DF$BlockOrder, "'1'='FirstCatBApp'; '2'='FirstCatAApp'") # Defining the "Subject", "Stimuli", "Movement" and "Category" as factors. DF$Subject <- factor(DF$pp) DF$Stimuli <- factor(DF$Stimuli) DF$Movement <- factor(DF$Movement) DF$Category <- factor(DF$Category) ### Cleaning dataset and Data exclusion ------------------------------------------------------------ # We subset the training phases DF <- subset (DF, BlockName=="blkTestCatAAp" | BlockName=="blkTestCatBAp") # We subset participants with screen resolution < 1200*675 DF <- subset (DF, screenWidth>1200) # OPTIONAL: We keep only participants being fluent in English #DF <- subset (DF, language==1) # Average response time ------------------------------------------------------ # Calculating the average response time per participant DFRTm <- data.frame (DF$RT, DF$Subject) DFRTm <- aggregate (DFRTm[, 1], list(DFRTm$DF.Subject), mean) DFRTm <- rename.vars (DFRTm, c("Group.1","x"), c("Subject","RT_Mean")) arrange(DFRTm,RT_Mean) # Subsetting participants having unusual mean response time DF <- merge (DF, DFRTm, by.x.=c(Subject), all=TRUE) DF<- subset(DF, RT_Mean < 3000) # Incorrect trials ------------------------------------------------------------ # Calculating the rate of incorrect trials per participant DFACC <- data.frame (DF$ACC, DF$Subject) DFACC <- aggregate (DFACC[, 1], list(DFACC$DF.Subject), sum) DFACC <- rename.vars (DFACC, c("Group.1","x"), c("Subject","ACC_Sum")) DFACC$ACC_rate <- (DFACC$ACC_Sum * 100) /80 arrange(DFACC,ACC_rate) # Rate of participants having less than 60% of correct trials # NB : we chose 60 % because under this limit, impossible to compute mean per condition (nrow(DFACC[DFACC$ACC_rate <60, ]) / nrow(DFACC) ) # Subsetting participants having less than 60% of correct trials DF <- merge (DF, DFACC, by.x.=c(Subject), all=TRUE) DF<- subset(DF, ACC_rate>60) # Overall rate of incorrect trials 1 - (nrow(DF[DF$ACC == 1, ]) / nrow(DF) ) # Subsetting incorrect trials. DF <- DF[DF$ACC == 1, ] # RT Cutoff ------------------------------------------------------------------- # Testing different cut-off to have the best RT distribution # define RT filters (minimum and maximum). #min <- 300 #max <- 1000 #min <- 350 #max <- 1500 #min <- 400 #max <- 2000 min <- 450 max <- 2500 #min <- 500 #max <- 3000 # Rate of exclusion if we keep only the RTs above min and below max 1 - (nrow(DF[DF$RT > min & DF$RT < max, ]) / nrow(DF) ) # We subset RTs above min and below max DF <- DF[DF$RT > min & DF$RT < max, ] # DV transformation ----------------------------------------------------------- # Inverse transformation of RTs DF$RTinv <- 1/DF$RT # Logarithmic transformation of RTs #DF$RTl <- log(DF$RT) ## Checking for the best filter according to the resulting distribution ----------------------------------------------------------- # With RT #ggplot(DF,aes(RT)) + # geom_histogram(binwidth=100) # With RTinv #ggplot(DF,aes(RTinv)) + # geom_histogram(binwidth=.0001) # With RTlog #ggplot(DF,aes(RTl)) + # geom_histogram(binwidth=.050) ### Recoding factors ------------------------------------------------------------ DF <- within(DF, { Category_C <- -0.5 * (Category == "catB") + 0.5 * (Category == "catA") Category_CatA <- 1 * (Category == "catB") + 0 * (Category == "catA") Category_CatB <- 0 * (Category == "catB") + 1 * (Category == "catA") Movement_C <- -0.5 * (Movement == "avoid") + 0.5 * (Movement == "appr") Movement_App <- 1 * (Movement == "avoid") + 0 * (Movement == "appr") Movement_Av <- 0 * (Movement == "avoid") + 1 * (Movement == "appr") BlockOrder_C <- -0.5 * (BlockOrder == "FirstCatAApp") + 0.5 * (BlockOrder == "FirstCatBApp") Gender_C <- -0.5 * (gender == 1) + 0.5 * (gender == 2) }) ### Mixed-modeling -------------------------------------------------------------- # Compatibility effect -------------------------------------------------------- # Loading the lmerTest package now (to avoid conflit with the dplyr package used before) library(lmerTest) # Movement by Category interaction ------------------------------------------------------ fitA.lmer <- lmer(RT ~ Movement_C * Category_C + (Movement_C * Category_C| Subject) + (Movement_C | Stimuli), data = DF) summary(fitA.lmer) # Block Order moderation ------------------------------------------------------ fitA.lmer <- lmer(RTinv ~ Movement_C * Category_C * BlockOrder_C + (Movement_C * Category_C | Subject) + (Movement_C | Stimuli), data = DF) summary(fitA.lmer) # Simple effect analysis ------------------------------------------------------ # Difference between approach and avoidance for stimuli of Category A ---------------- fitA.lmer <- lmer(RTinv ~ Movement_C * Category_CatA + (Movement_C * Category_CatA| Subject) + (Movement_C | Stimuli) , data = DF) summary(fitA.lmer) # Difference between approach and avoidance for stimuli of Category B ---------------- fitA.lmer <- lmer(RTinv ~ Movement_C * Category_CatB + (Movement_C * Category_CatB| Subject) + (Movement_C | Stimuli) , data = DF) summary(fitA.lmer) # Difference between Category A and Category B for approach ------------------------- fitA.lmer <- lmer(RTinv ~ Movement_App * Category_C + (Movement_App * Category_C| Subject) + (Movement_App | Stimuli) , data = DF) summary(fitA.lmer) # Difference between Category A and Category B for avoidance ------------------------ fitA.lmer <- lmer(RTinv ~ Movement_Av * Category_C + (Movement_Av * Category_C| Subject) + (Movement_Av | Stimuli) , data = DF) summary(fitA.lmer) # Descriptive statistics ------------------------------------------------------ # Data wrangling -------------------------------------------------------------- DF_short <- cast(DF, Subject ~ Movement + Category, value = "RT", # to be rerun with "RT" to compute the means and se mean, fill = NA) # New variables --------------------------------------------------------------- DF_short <- within(DF_short, { Comp <- (appr_catA + avoid_catB) / 2 Incomp <- (appr_catB + avoid_catA) / 2 diff <- Incomp - Comp MeanAppr <- (appr_catA + appr_catB) / 2 MeanAvoid <- (avoid_catA + avoid_catB) / 2 MeanCatA <- (appr_catA + avoid_catA) / 2 MeanCatB <- (appr_catB + avoid_catB) / 2 diffMov <- MeanAppr - MeanAvoid diffVal <- MeanCatA - MeanCatB diffApVal <- appr_catA - appr_catB diffAvVal <- avoid_catA - avoid_catB diffMovCatA <- appr_catA - avoid_catA diffMovCatB <- appr_catB - avoid_catB Pourc <- ifelse(diff < 0 , "N", "O") }) # Means per condition. attach(DF_short) MeanCompatibility <- colMeans(cbind(Incomp, Comp), na.rm=TRUE) MeansCondition <- colMeans(cbind(appr_catA, appr_catB, avoid_catA, avoid_catB),na.rm=TRUE ) ci.fun(appr_catA) ci.fun(appr_catB) ci.fun(avoid_catA) ci.fun(avoid_catB) detach(DF_short) # SE per condition. sd(DF_short$MeanAppr) / sqrt(length(DF_short$MeanAppr)) sd(DF_short$MeanAvoid) / sqrt(length(DF_short$MeanAvoid)) sd(DF_short$MeanCatA) / sqrt(length(DF_short$MeanCatA)) sd(DF_short$MeanCatB) / sqrt(length(DF_short$MeanCatB)) sd(DF_short$appr_catA) / sqrt(length(DF_short$appr_catA)) sd(DF_short$appr_catB) / sqrt(length(DF_short$appr_catB)) sd(DF_short$avoid_catA) / sqrt(length(DF_short$avoid_catA)) sd(DF_short$avoid_catB) / sqrt(length(DF_short$avoid_catB)) # dz effect size Category*Movement --------------------------------------------------------------------- dz <- mean(DF_short$diff) / sd(DF_short$diff) # dz effect size Movement --------------------------------------------------------------------- dz <- mean(DF_short$diffMov) / sd(DF_short$diffMov) # dz effect size Category --------------------------------------------------------------------- dz <- mean(DF_short$diffVal) / sd(DF_short$diffVal) # dz effect size Approach_Category --------------------------------------------------------------------- dz <- mean(DF_short$diffApVal) / sd(DF_short$diffApVal) # dz effect size Avoid_Category --------------------------------------------------------------------- dz <- mean(DF_short$diffAvVal) / sd(DF_short$diffAvVal) # dz effect size Movement_CatA --------------------------------------------------------------------- dz <- mean(DF_short$diffMovCatA) / sd(DF_short$diffMovCatA) # dz effect size Movement_CatB --------------------------------------------------------------------- dz <- mean(DF_short$diffMovCatB) / sd(DF_short$diffMovCatB) # Creating a figure --------------------------------------------------------------------- # Creating a data table with the results summaries (mean, CI) # Replace means and CIs in blue by means and CIs computed above DFf2<-tribble( ~Movement, ~Category, ~mean, ~CI, "Approach","CategoryA", 892.7671, 53.0205, "Approach","CategoryB", 973.1305, 52.8475, "Avoidance","CategoryA", 956.3348, 58.7052, "Avoidance","CategoryB", 905.8374, 47.2511 ) # From the data table, we create a figure ggplot(DFf2, aes(x = Movement,y = mean,fill=Category)) + geom_bar(position="dodge",stat="identity",width=0.6) + scale_x_discrete(limits=c("Approach","Avoidance")) + theme_bw() + theme(panel.border = element_blank(), panel.grid.major.x = element_blank(), axis.text.x = element_text(size=16,color = "black"), axis.title.y = element_text(size=16), axis.ticks.y = element_blank(), axis.text.y = element_text(size=16,color = "black"), axis.ticks.x = element_blank(), axis.line = element_line(colour="black",size=0.2), legend.text = element_text(size=16)) + coord_cartesian(ylim=c(800,1020)) + geom_errorbar(aes(ymin=mean-CI,ymax=mean+CI), position=position_dodge(0.6),width=0.1) + labs(x="",y="Response time (ms)") + guides(fill=guide_legend((title=NULL))) + scale_fill_manual(values=c("grey70","grey40")) # To save the figure in the R repertoire ggsave("Figure.tiff", dpi=80, width = 20, height = 20, units = "cm") # Demographic informations ---------------------------------------------------- # Data wrangling -------------------------------------------------------------- DF_demo <- cast(DF, pp + age + gender ~ Movement, value = "RT", mean, fill = NA) # Age ------------------------------------------------------------------------- mean(DF_demo$age) sd(DF_demo$age) # Gender ---------------------------------------------------------------------- # (1 = male) table(DF_demo$gender)