MatchIdx = match(SubjSess,CondData$PARTICIPANTS_SESSIONS)
CountRVIPRate$Conditions = CondData$CONDITIONS[MatchIdx]
# Load all data (for each trial)
RVIPDataAll = readRDS(paste(RelativePath,"Data/RVIPDataAll.rds",sep="/"))
RVIPDataAll$Participant = as.integer(str_replace_all(RVIPDataAll$Participant,"P",""))
RVIPDataAll$Participant = as.factor(RVIPDataAll$Participant)
RVIPDataAll$Session = as.factor(RVIPDataAll$Session)
# Sort by natural order
RVIPDataAll = RVIPDataAll[stri_order(RVIPDataAll$Participant, numeric = TRUE),]
# Then remove trial outliers blocks for each subject/session
BlocOutliers = which(RVIPDataRate$Response_thr == "FA.Outlier" | RVIPDataRate$Response_thr == "MISS.Outlier")
Idx = NA
for (k in 1:length(BlocOutliers)) {
PartNum = str_replace(RVIPDataRate$Participant[BlocOutliers[k]],"P","")
SessionNum = RVIPDataRate$Session[BlocOutliers[k]]
BlocNum = RVIPDataRate$bloc[BlocOutliers[k]]
# Idx to remove the outliers
# Idx to remove the outliers
Idx = c(Idx, which(RVIPDataAll$Participant==PartNum & RVIPDataAll$Session == SessionNum & RVIPDataAll$bloc == BlocNum))
}
Idx = unique(Idx[-1])
# Remove the outliers
RVIPDataAll = RVIPDataAll[-Idx,]
# Add the conditions column
RVIPDataAll$Conditions = rep(NA,dim(RVIPDataAll)[1])
SubjSess = paste("P",RVIPDataAll$Participant,"_",RVIPDataAll$Session, sep="")
for (k in 1:length(CondData$PARTICIPANTS_SESSIONS)) {
Idx = which(SubjSess == CondData$PARTICIPANTS_SESSIONS[k])
RVIPDataAll$Conditions[Idx] = rep(CondData$CONDITIONS[k],length(Idx))
}
# REJECTING OUTLIER TRIALS AND RECOMPUTING STATISTICAL MEASURES
# Count number of rejected participants
SummaryRVIP = as.integer(str_remove_all(CountRVIPRate$Participant[CountRVIPRate$Kept<=40],"P"))
SummaryRVIPCond = CountRVIPRate$Conditions[CountRVIPRate$Kept<=40]
# Replacing rejected participants data with NAs
# https://stackoverflow.com/questions/16844613/aggregate-methods-treat-missing-values-na-differentl
Idx = NA
if (length(SummaryRVIP)>0){
for (k in 1:length(SummaryRVIP)){ # If participants/sessions need to be ignored
Idx[k] = which(RVIPData$Participant == SummaryRVIP[k] & RVIPData$Conditions == SummaryRVIPCond[k])
}
} else {
Idx = seq(1,dim(GNGData)[1])*-1 # This will cancel the effect of the -Idx below
}
RVIPData[Idx,seq(3,dim(RVIPData)[2]-1)] = NA
# Recompute the statistical measures and replace in original data frame (i.e. RVIPData)
# meanRT_clean
TEMPDat = aggregate(RVIPDataAll$RT_HIT_Clean, list(RVIPDataAll$Participant, RVIPDataAll$Session), mean, na.rm = T, na.action="na.pass")
TEMPDat = TEMPDat[stri_order(TEMPDat$Group.1, numeric = TRUE),] # Sort natural
RVIPData$meanRT_clean[-Idx] = TEMPDat[-Idx,dim(TEMPDat)[2]]
# medianRT_clean
TEMPDat = aggregate(RVIPDataAll$RT_HIT_Clean, list(RVIPDataAll$Participant, RVIPDataAll$Session), median, na.rm = T, na.action="na.pass")
TEMPDat = TEMPDat[stri_order(TEMPDat$Group.1, numeric = TRUE),] # Sort natural
RVIPData$medianRT_clean[-Idx] = TEMPDat[-Idx,dim(TEMPDat)[2]]
# stdRT_clean
TEMPDat = aggregate(RVIPDataAll$RT_HIT_Clean, list(RVIPDataAll$Participant, RVIPDataAll$Session), sd, na.rm = T)
TEMPDat = TEMPDat[stri_order(TEMPDat$Group.1, numeric = TRUE),] # Sort natural
RVIPData$stdRT_clean[-Idx] = TEMPDat[-Idx,dim(TEMPDat)[2]]
# countHIT
TEMPDat = aggregate(RVIPDataAll$HIT, list(RVIPDataAll$Participant, RVIPDataAll$Session), sum, na.rm = T)
TEMPDat = TEMPDat[stri_order(TEMPDat$Group.1, numeric = TRUE),] # Sort natural
RVIPData$countHIT[-Idx] = TEMPDat[-Idx,dim(TEMPDat)[2]]
# countMISS
TEMPDat = aggregate(RVIPDataAll$MISS, list(RVIPDataAll$Participant, RVIPDataAll$Session), sum, na.rm = T)
TEMPDat = TEMPDat[stri_order(TEMPDat$Group.1, numeric = TRUE),] # Sort natural
RVIPData$countMISS[-Idx] = TEMPDat[-Idx,dim(TEMPDat)[2]]
# countFA
TEMPDat = aggregate(RVIPDataAll$FA, list(RVIPDataAll$Participant, RVIPDataAll$Session), sum, na.rm = T)
TEMPDat = TEMPDat[stri_order(TEMPDat$Group.1, numeric = TRUE),] # Sort natural
RVIPData$countFA[-Idx] = TEMPDat[-Idx,dim(TEMPDat)[2]]
# %HIT
RVIPData$`%HIT`[-Idx] = (RVIPData$countHIT/(RVIPData$countHIT+RVIPData$countMISS)*100)[-Idx]
# %MISS
RVIPData$`%MISS`[-Idx] = (RVIPData$countMISS/(RVIPData$countHIT+RVIPData$countMISS)*100)[-Idx]
# %FA
RVIPData$`%FA`[-Idx]= (RVIPData$countFA/192*100)[-Idx]
# A. Questionnaires data
QuestDataTemp = Data[,c(1,seq(48,61),63)]
QuestData = QuestDataTemp[Data$Session=="S1",-c(1,dim(QuestDataTemp)[2]-2)] # This will be the final dataset
# Averaging over the 3 measures of STAI_TRAIT
STAITrait = QuestDataTemp[,c(1,15)] %>%  group_by(Subjects) %>% summarise_each(funs( mean( .,na.rm = TRUE)))
QuestDataTemp = QuestDataTemp[c(1,dim(QuestDataTemp)[2],dim(QuestDataTemp)[2]-2)]
# # need to rearrange the columns so that STAI_STATE is divided into 3 columns (for the 3 sessions)
STAIState = QuestDataTemp %>% spread(Conditions, STAI_State_Global)
colnames(STAIState)[c(2,3,4)] = paste("STAI_State",sep="_",colnames(STAIState)[c(2,3,4)])
QuestData = cbind(QuestData,STAIState[,c(2,3,4)])
QuestData$STAI_Trait_Global = STAITrait$STAI_Trait_Global
QuestData$Subjects = unique(Data$Subjects)
QuestData = QuestData[,-which(colnames(QuestData)=="Conditions")] # Removing the conditions column
# Computing MAD
QuestOutliers = ComputeMAD(QuestData)
# Winsorization
WinsorQuestData = QuestData[,-which(colnames(QuestData)=="Subjects")]
for (k in 1:length(QuestOutliers[[2]])) {
if (is.integer(WinsorQuestData[,k])==T) {
# Round if data type is integer
Bounds = as.integer(round(QuestOutliers[[2]][[k]]))
}  else {
Bounds = QuestOutliers[[2]][[k]]
}
Bounds = ifelse(Bounds < 0, 0, Bounds) # negative values wouldn't make any sense
# Exclude the case where bounds are [0 0]
if (length(which( Bounds!=0)) >= 1){
WinsorQuestData[,k] = Winsorize(QuestData[,k], minval = Bounds[1],maxval = Bounds[2], probs = c(0.05, 0.95))
}
}
WinsorQuestData$Subjects = QuestData$Subjects
# Gathering data to go back to having a "Condition" columns
WinsorQuestDataStats = Data;
VarList = colnames(WinsorQuestData)[-seq(which(colnames(WinsorQuestData)=="STAI_State_GC-TD"),dim(WinsorQuestData)[2])]
for (k in 1:length(VarList)) {
for (j in 1:dim(WinsorQuestDataStats)[1]) {
WinsorQuestDataStats[WinsorQuestDataStats$Session=="S1",which(colnames(WinsorQuestDataStats)==VarList[k])] = WinsorQuestData[,which(colnames(WinsorQuestData)==VarList[k])]
}
}
# Only for STAI
Cols = c("STAI_State_GC-TD","STAI_State_GD-TC","STAI_State_GD-TD")
WinsorQuestData.long = WinsorQuestData[,c(17,seq(14,16))] %>% gather(Conditions, STAI_State_Global, all_of(Cols))
SortIdx = order(WinsorQuestData.long$Subjects)
WinsorQuestData.long = WinsorQuestData.long[SortIdx,]
WinsorQuestData.long$Conditions = as.factor(WinsorQuestData.long$Conditions)
levels(WinsorQuestData.long$Conditions) = c(CondAssign[1],CondAssign[2],CondAssign[3])
for (j in 1:dim(WinsorQuestDataStats)[1]) {
# Find position
Idx = which(WinsorQuestData.long$Subjects == WinsorQuestDataStats$Subjects[j] &
WinsorQuestData.long$Conditions ==  WinsorQuestDataStats$Conditions[j])
if (length(Idx)==1) {
# Replace in original data frame
WinsorQuestDataStats[j,which(colnames(WinsorQuestDataStats)=="STAI_State_Global")] = WinsorQuestData.long[Idx,dim(WinsorQuestData.long)[2]]
}
}
# Saving as .xlsx in two separate sheets (outlier data and winsorized data)
FilePath = paste(RelativePath,"WinsorData", paste("CAFQuestionnaires_Full_",DateTime,".xlsx",sep=""),sep="/")
ListDFs = list("Outliers" = Data, "Winsorized" = WinsorQuestDataStats, "Final" = Data) # The final sheet can be manually adjusted before analyses
write.xlsx(ListDFs, file = FilePath)
# B. Cognitive tasks behavioral data
CogData = data.frame(Subjects = BehavData$Subjects[BehavData$GNG_HITMISS=="HIT"],
Conditions = BehavData$Conditions[BehavData$GNG_HITMISS=="HIT"],
RVIPMeanRT = RVIPData$meanRT_clean,
RVIPHIT = RVIPData$`%HIT`,
RVIPFA = RVIPData$`%FA`,
RVIPMISS = RVIPData$`%MISS`,
GNGMeanRT = GNGData$meanRT_clean,
GNGHIT = GNGData$`%HIT`,
GNGCR = GNGData$`%CR`,
GNGFA = GNGData$`%FA`,
GNGMISS = GNGData$`%MISS`)
# Reorder data frame according to Subjects and Conditions
CogData = CogData[order(CogData[,1],CogData[,2]),]
NEWCogData = data.frame(Subjects = unique(CogData$Subjects))
VarList = colnames(CogData)[-c(1,2)]
for (k in 1:length(VarList)) {
Temp = CogData[,c(1,2,which(colnames(CogData) == VarList[k]))] %>% spread(Conditions, VarList[k])
colnames(Temp)[c(2,3,4)] = paste(VarList[k],colnames(Temp)[c(2,3,4)],sep="_")
NEWCogData = cbind(NEWCogData,Temp[,c(2,3,4)])
}
BehavOutliers = ComputeMAD(NEWCogData)
# Winsorization
NEWCogDataNoSubj = NEWCogData[,-which(colnames(NEWCogData)=="Subjects")]
WinsorCogData = NEWCogDataNoSubj
for (k in 1:length(BehavOutliers[[2]])) {
if (is.integer(WinsorCogData[,k])==T) {
# Round if data type is integer
Bounds = as.integer(round(BehavOutliers[[2]][[k]]))
}  else {
Bounds = BehavOutliers[[2]][[k]]
}
Bounds = ifelse(Bounds < 0, 0, Bounds) # negative values wouldn't make any sense
# Exclude the case where bounds are [0 0]
if (length(which( Bounds!=0)) >= 1){
WinsorCogData[,k] = Winsorize(NEWCogDataNoSubj[,k], minval = Bounds[1],maxval = Bounds[2], probs = c(0.05, 0.95))
}
}
WinsorCogData$Subjects = NEWCogData$Subjects
# Gathering data to go back to having a "Condition" columns
WinsorCogDataStats = CogData
for (k in 1:length(VarList)) {
VarToGather = paste(VarList[k],sort(unique(CogData$Conditions)),sep="_")
Idx = which(colnames(WinsorCogData)==VarToGather[1])
WinsorCogData.long = WinsorCogData[,c(Idx,Idx+1,Idx+2,28)] %>% gather(Conditions, !!VarList[k], all_of(VarToGather))
# !! see: https://stackoverflow.com/questions/57677444/use-of-gather-in-a-loop-key-and-value-as-a-variable
SortIdx = order(WinsorCogData.long$Subjects)
WinsorCogData.long = WinsorCogData.long[SortIdx,]
# replace Variables string in Conditions
WinsorCogData.long$Conditions = str_replace_all(WinsorCogData.long$Conditions,paste(VarList[k],"_",sep=""),"")
# Reorder data frame according to Subjects and Conditions
WinsorCogData.long = WinsorCogData.long[order(WinsorCogData.long[,1],WinsorCogData.long[,2]),]
# Replace in original data frame
WinsorCogDataStats[,which(colnames(WinsorCogDataStats)==VarList[k])] = WinsorCogData.long[,dim(WinsorCogData.long)[2]]
}
apply(WinsorCogDataStats, 2, function(x) {ifelse(x < 0, 0, x)})
# Saving as .xlsx in two separate sheets (outlier data and winsorized data)
FilePath = paste(RelativePath,"WinsorData", paste("CAFCognitive_Full_",DateTime,".xlsx",sep=""),sep="/")
ListDFs = list("Outliers" = CogData, "Winsorized" = WinsorCogDataStats, "Final" = CogData) # The final sheet can be manually adjusted before analyses)
write.xlsx(ListDFs, file = FilePath)
# Output files
OutFiles = c("CAFEEG_GNG_N2_Full_","CAFEEG_GNG_P3_Full_","CAFEEG_RVIP_P3_Full_")
for (j in 1:3) {
# Reorder data frame according to Subjects and Levels
EEGDataList[[j]] = EEGDataList[[j]][order(EEGDataList[[j]][,9],EEGDataList[[j]][,2]),]
NEWEEGData = data.frame(Subjects = unique(EEGDataList[[j]]$Subjects))
VarList = colnames(EEGDataList[[j]])[-c(1,2,9)]
for (k in 1:length(VarList)) {
Temp = EEGDataList[[j]][,c(9,2,which(colnames(EEGDataList[[j]]) == VarList[k]))] %>% spread(Levels, VarList[k])
colnames(Temp)[c(2,3,4)] = paste(VarList[k],colnames(Temp)[c(2,3,4)],sep="_")
NEWEEGData = cbind(NEWEEGData,Temp[,c(2,3,4)])
}
BehavOutliers = ComputeMAD(NEWEEGData)
# Winsorization
WinsorEEGData = NEWEEGData[,-1]
for (k in 1:length(BehavOutliers[[2]])) {
if (is.integer(WinsorEEGData[,k])==T) {
# Round if data type is integer
Bounds = as.integer(round(BehavOutliers[[2]][[k]]))
}  else {
Bounds = BehavOutliers[[2]][[k]]
}
# Exclude the case where bounds are [0 0]
if (length(which( Bounds!=0)) >= 1){
WinsorEEGData[,k] = Winsorize(WinsorEEGData[,k], minval = Bounds[1],maxval = Bounds[2], probs = c(0.05, 0.95))
}
}
WinsorEEGData$Subjects = NEWEEGData$Subjects
# Gathering data to go back to having a "Condition" columns
WinsorEEGDataStats = EEGDataList[[j]]
for (k in 1:length(VarList)) {
VarToGather = paste(VarList[k],sort(unique(EEGDataList[[j]]$Levels)),sep="_")
Idx = which(colnames(WinsorEEGData)==VarToGather[1])
WinsorEEGData.long = WinsorEEGData[,c(Idx,Idx+1,Idx+2,19)] %>% gather(Conditions, !!VarList[k], all_of(VarToGather))
# !! see: https://stackoverflow.com/questions/57677444/use-of-gather-in-a-loop-key-and-value-as-a-variable
SortIdx = order(WinsorEEGData.long$Subjects)
WinsorEEGData.long = WinsorEEGData.long[SortIdx,]
# replace Variables string in Conditions
WinsorEEGData.long$Conditions = str_replace_all(WinsorEEGData.long$Conditions,paste(VarList[k],"_",sep=""),"")
# Reorder data frame according to Subjects and Conditions
WinsorEEGData.long = WinsorEEGData.long[order(WinsorEEGData.long[,1],WinsorEEGData.long[,2]),]
# Replace in original data frame
WinsorEEGDataStats[,which(colnames(WinsorEEGDataStats)==VarList[k])] = WinsorEEGData.long[,dim(WinsorEEGData.long)[2]]
}
# Saving as .xlsx in two separate sheets (outlier data and winsorized data)
FilePath = paste(RelativePath,"WinsorData", paste(OutFiles[j],DateTime,".xlsx",sep=""),sep="/")
ListDFs = list("Outliers" = EEGDataList[[j]], "Winsorized" = WinsorEEGDataStats, "Final" = EEGDataList[[j]]) # The final sheet can be manually adjusted before analyses)
write.xlsx(ListDFs, file = FilePath)
}
# Separated by conditions
Plots = list(); Conditions = c(CondAssign[1],CondAssign[2],CondAssign[3]);
for (k in seq(1,3)) {
TEMPGNGDataRate = GNGDataRate[GNGDataRate$Conditions==Conditions[k],] # One condition at a time
TEMPGNGDataRate = TEMPGNGDataRate[TEMPGNGDataRate$Response=="FA" | TEMPGNGDataRate$Response=="MISS",] # Only keeping lines related to FA / MISS
Plots[[k]] = MainPlot(TEMPGNGDataRate,TEMPGNGDataRate$bloc,TEMPGNGDataRate$rate,TEMPGNGDataRate$Response,TEMPGNGDataRate$logical_thr,paste("GNG response rate:",Conditions[k]))
plot(Plots[[k]])
}
# Sorting according to natural order
PlotBar(CountGNGRate,CountGNGRate$Participant,CountGNGRate$Kept,"GNG kept blocs (in %)",40,T)
# Separated by conditions
Plots = list(); Conditions = c(CondAssign[1],CondAssign[2],CondAssign[3]);
for (k in seq(1,3)) {
TEMPRVIPDataRate = RVIPDataRate[RVIPDataRate$Conditions==Conditions[k],] # One condition at a time
TEMPRVIPDataRate = TEMPRVIPDataRate[TEMPRVIPDataRate$Response=="FA" | TEMPRVIPDataRate$Response=="MISS",] # Only keeping lines related to FA / MISS
Plots[[k]] = MainPlot(TEMPRVIPDataRate,TEMPRVIPDataRate$bloc,TEMPRVIPDataRate$rate,TEMPRVIPDataRate$Response,TEMPRVIPDataRate$logical_thr,paste("RVIP response rate:",Conditions[k]))
plot(Plots[[k]])
}
# Sorting according to natural order
PlotBar(CountRVIPRate,CountRVIPRate$Participant,CountRVIPRate$Kept,"RVIP kept blocs (in %)",40,T)
# Extract the replies to the open question
GNGReply =readRDS(paste(RelativePath,"Data/GNGReply.rds",sep="/"))
# Finding matching patterns
GNGPatternsList = FindMatch(GNGReply, c(""), "GNG", CondData)
GNGOpenQuestPercDF = GNGPatternsList[[length(GNGPatternsList)-3]]
GNGOpenQuestBoulDF= GNGPatternsList[[length(GNGPatternsList)-2]]
plot(GNGPatternsList[[length(GNGPatternsList)]])
# Table with participants' response
GNGReplyDF = data.frame(SbjSession = names(GNGReply), Resp = rep(NA,length(GNGReply)))
GNGReplyDF[,2] = data.frame(matrix(unlist(GNGReply), nrow=length(GNGReply), byrow=T))
IdxNatSort = str_order(GNGReplyDF$SbjSession, numeric = TRUE)
GNGReplyDF = GNGReplyDF[IdxNatSort,]
# Print table
datatable(GNGReplyDF, options = list(rows.print = 15), caption = htmltools::tags$caption(
style = 'font-size:150%',"GNG: Participants' responses to open questions"))
# Extract the replies to the open question
RVIPReply =readRDS(paste(RelativePath,"Data/RVIPReply.rds",sep="/"))
# Finding matching patterns
RVIPPatternsList = FindMatch(RVIPReply, c(""), "RVIP", CondData)
RVIPOpenQuestPercDF = RVIPPatternsList[[length(RVIPPatternsList)-3]]
RVIPOpenQuestBoulDF= RVIPPatternsList[[length(RVIPPatternsList)-2]]
plot(RVIPPatternsList[[length(RVIPPatternsList)]])
# Table with participants' response
RVIPReplyDF = data.frame(SbjSession = names(RVIPReply), Resp = rep(NA,length(RVIPReply)))
RVIPReplyDF[,2] <- data.frame(matrix(unlist(RVIPReply), nrow=length(RVIPReply), byrow=T))
IdxNatSort = str_order(RVIPReplyDF$SbjSession, numeric = TRUE)
RVIPReplyDF = RVIPReplyDF[IdxNatSort,]
# Print table
datatable(RVIPReplyDF, options = list(rows.print = 15), caption = htmltools::tags$caption(
style = 'font-size:150%',"RVIP: Participants' responses to open questions"))
# Global score for items 1+2+6
DPEQCriteria = Data$DPEQ1_Goal + Data$DPEQ2_Doubts + Data$DPEQ6_Placebo
DPEQData = cbind(Data,DPEQCriteria)
PlotBar(DPEQData,DPEQData$Subjects,DPEQData$DPEQCriteria,"DPEQ score (items 1+2+6)",18,F)
# Separated questions
DPEQData2 = data.frame (Data = c(Data$DPEQ1_Goal,Data$DPEQ2_Doubts,Data$DPEQ6_Placebo),
Session = rep(c("Item1 Goal","Item2 Doubts","Item6 Placebo"),each = dim(Data)[1]),
Subjects = rep(Data$Subjects,3))
DPEQData2 = DPEQData2[!(is.na(DPEQData2[,1])),] # Remove NA lines
DPEQData2$Session = as.factor(DPEQData2$Session)
ggplot(DPEQData2,aes(x=Subjects, y=Data, fill=Subjects, label=round(Data, digits = 2))) +
geom_bar(stat = "identity") +
labs(x = "Subjects", y = "Scores") +
theme_bw() +
theme(legend.position = "none",
axis.title.y = element_text(size=12*TxtS),
axis.title.x = element_text(size=12*TxtS),
axis.text.x = element_text(size=10*TxtS, angle = 90),
axis.text.y = element_text(size=10*TxtS),
legend.text = element_text(size=10*TxtS),
strip.text.x = element_text(size = 12*TxtS),
plot.title = element_text(size=14*TxtS,hjust = 0.5)) + facet_grid(.~DPEQData2$Session)
# Print table
TEMPDF = Data[seq(3,dim(Data)[1],3),c(1,25,27)]
colnames(TEMPDF) = c("Subjects","Pensez-vous avoir deviné le but ?","Avez-vous eu des doutes ? ")
datatable(TEMPDF, options = list(rows.print = 15), caption = htmltools::tags$caption(
style = 'font-size:150%',"DPEQ debriefing questionnaire"))
# Plot the number of retained trials after pre-processing
ggplot(ERPData, aes(x = Subj, y = N_Trials_in_ERP, fill = Trigger)) +
geom_bar(stat = "identity", position=position_dodge()) +
labs(title = "Results of ERP pre-processing",
subtitle = paste("The horizontal dashed line represent the limit of",ERPThresh,"trials"),
x = "Participant & Session", y = "Number of kept trials", fill = "Triggers") +
facet_grid(Trigger~Conditions, scales="free") + # .~Trigger # , ncol  = ceiling(NPlots/4)
geom_hline(yintercept=ERPThresh,linetype=8, color="black") +
theme_bw() +
theme(axis.text.x = element_text(size=10*TxtS, angle=90),
axis.text.y = element_text(size=10*TxtS),
axis.title.y = element_text(size=14*TxtS),
axis.title.x = element_text(size=14*TxtS),
strip.text.x = element_text(size=14*TxtS),
strip.text.y = element_text(size=14*TxtS),
plot.title = element_text(size=20*TxtS,hjust = 0.5),
plot.subtitle = element_text(size=16*TxtS,hjust = 0.5),
legend.title = element_text(size=14*TxtS),
legend.text = element_text(size=12*TxtS))
# Print table
datatable(ERPData, options = list(rows.print = 15), caption = htmltools::tags$caption(
style = 'font-size:150%','Results of ERP pre-processing (autoERP)'))
# Find list of participant with ERP trials < ERPThresh
LowTrialsIdx = ERPData$N_Trials_in_ERP<ERPThresh
LowFileNames = ERPData$File_Name[LowTrialsIdx]
LowTriggers = ERPData$Trigger[LowTrialsIdx]
LowConditions = ERPData$Conditions[LowTrialsIdx]
# Manually adding P900 missing data for GNG-CR & GNG-HIT in GD-TD condition
LowFileNames = c(LowFileNames,"\\GNG\\CR\\P900_S2_GNG","\\GNG\\HIT\\P900_S2_GNG")
LowConditions = c(LowConditions,"GD-TD","GD-TD")
Levs = levels(LowTriggers)
LowTriggers = c(as.character(LowTriggers), "GNG_CR","GNG_HIT")
LowTriggers = factor(LowTriggers, levels = Levs)
Idx = str_order(LowFileNames, numeric = TRUE)
LowFileNames = LowFileNames[Idx]
LowTriggers = LowTriggers[Idx]
LowConditions = LowConditions[Idx]
# Building data frame
ERPRejSubj = data.frame(Names=rep(NA,length(LowFileNames)),
Triggers=rep(NA,length(LowTriggers)),
Conditions=LowConditions)
for (k in 1:length(LowFileNames)) {
ERPRejSubj$Names[k] = parse_number(LowFileNames[k])[1]
ERPRejSubj$Triggers[k] = as.character(LowTriggers[k])
}
# Summary of the rejected subjects for each condition and variable of interest
BlocGNG_C1 = as.integer(str_remove_all(CountGNGRate[CountGNGRate$Kept<=40 & CountGNGRate$Conditions=="GD-TD",1],"P"))
BlocGNG_C2 = as.integer(str_remove_all(CountGNGRate[CountGNGRate$Kept<=40 & CountGNGRate$Conditions=="GC-TD",1],"P"))
BlocGNG_C3 = as.integer(str_remove_all(CountGNGRate[CountGNGRate$Kept<=40 & CountGNGRate$Conditions=="GD-TC",1],"P"))
BlocRVIP_C1 = as.integer(str_remove_all(CountRVIPRate[CountRVIPRate$Kept<=40 & CountRVIPRate$Conditions=="GD-TD",1],"P"))
BlocRVIP_C2 = as.integer(str_remove_all(CountRVIPRate[CountRVIPRate$Kept<=40 & CountRVIPRate$Conditions=="GC-TD",1],"P"))
BlocRVIP_C3 = as.integer(str_remove_all(CountRVIPRate[CountRVIPRate$Kept<=40 & CountRVIPRate$Conditions=="GD-TC",1],"P"))
DPEQ = as.integer(as.character(sort(DPEQData$Subjects[which(DPEQData$DPEQCriteria>18)])))
ERPRejSubj_C1 = subset(ERPRejSubj, Conditions=="GD-TD")
GNG_CR_C1 = ERPRejSubj_C1$Names[ERPRejSubj_C1$Triggers=="GNG_CR"]
RVIP_HITN1_C1 = ERPRejSubj_C1$Names[ERPRejSubj_C1$Triggers=="RVIP_HITN1"]
ERPRejSubj_C2 = subset(ERPRejSubj, Conditions=="GC-TD")
GNG_CR_C2 = ERPRejSubj_C2$Names[ERPRejSubj_C2$Triggers=="GNG_CR"]
RVIP_HITN1_C2 = ERPRejSubj_C2$Names[ERPRejSubj_C2$Triggers=="RVIP_HITN1"]
ERPRejSubj_C3 = subset(ERPRejSubj, Conditions=="GD-TC")
GNG_CR_C3 = ERPRejSubj_C3$Names[ERPRejSubj_C3$Triggers=="GNG_CR"]
RVIP_HITN1_C3 = ERPRejSubj_C3$Names[ERPRejSubj_C3$Triggers=="RVIP_HITN1"]
Size = max(length(DPEQ),length(BlocGNG_C1),length(BlocGNG_C2),length(BlocGNG_C3),length(BlocRVIP_C1),
length(BlocRVIP_C2),length(BlocRVIP_C3),length(GNG_CR_C1),length(RVIP_HITN1_C1), length(GNG_CR_C2),
length(RVIP_HITN1_C2), length(GNG_CR_C3),length(RVIP_HITN1_C3))
ValidPart =  data.frame(DPEQ = rep(NA,Size), BlocGNG_C1 = rep(NA,Size), BlocGNG_C2 = rep(NA,Size), BlocGNG_C3 = rep(NA,Size),
BlocRVIP_C1 = rep(NA,Size), BlocRVIP_C2 = rep(NA,Size), BlocRVIP_C3 = rep(NA,Size),
GNG_CR_C1 = rep(NA,Size), GNG_CR_C2 = rep(NA,Size), GNG_CR_C3 = rep(NA,Size),
RVIP_HITN1_C1 = rep(NA,Size), RVIP_HITN1_C2 = rep(NA,Size), RVIP_HITN1_C3 = rep(NA,Size))
for (k in 1:Size) {
ValidPart$DPEQ[k] = DPEQ[k]
ValidPart$BlocGNG_C1[k] = BlocGNG_C1[k]
ValidPart$BlocGNG_C2[k] = BlocGNG_C2[k]
ValidPart$BlocGNG_C3[k] = BlocGNG_C3[k]
ValidPart$BlocRVIP_C1[k] = BlocRVIP_C1[k]
ValidPart$BlocRVIP_C2[k] = BlocRVIP_C2[k]
ValidPart$BlocRVIP_C3[k] = BlocRVIP_C3[k]
ValidPart$GNG_CR_C1[k] = GNG_CR_C1[k]
ValidPart$GNG_CR_C2[k] = GNG_CR_C2[k]
ValidPart$GNG_CR_C3[k] = GNG_CR_C3[k]
ValidPart$RVIP_HITN1_C1[k] = RVIP_HITN1_C1[k]
ValidPart$RVIP_HITN1_C2[k] = RVIP_HITN1_C2[k]
ValidPart$RVIP_HITN1_C3[k] = RVIP_HITN1_C3[k]
}
# Counting the number of subjects which can be included for each condition and variable of interest
CompleteDPEQ = sum(DPEQData$DPEQCriteria<=18,na.rm=T)
CompleteBehavGNG = CountGNGRate %>% group_by(Conditions) %>% dplyr::count(Kept>40)
CompleteBehavGNG = CompleteBehavGNG[CompleteBehavGNG$`Kept > 40`==TRUE,]
CompleteBehavRVIP = CountRVIPRate %>% group_by(Conditions) %>% dplyr::count(Kept>40)
CompleteBehavRVIP = CompleteBehavRVIP[CompleteBehavRVIP$`Kept > 40`==TRUE,]
CompleteERP = ERPData %>% group_by(Trigger,Conditions) %>% dplyr::count(N_Trials_in_ERP>=ERPThresh)
CompleteERP = CompleteERP[CompleteERP$`N_Trials_in_ERP >= ERPThresh`==TRUE,]
# Already compensating for missing GNG for P900 in GD-TD condition since 43 subj in total in this condition (while 44 in the 2 others)
# Combining all together
CondOrder = c(3,1,2) # Reorder according the conditions order (since tibble will order alphabetically)
CompleteData = c(CompleteDPEQ, CompleteBehavGNG$n[CondOrder], CompleteBehavRVIP$n[CondOrder],
CompleteERP$n[CondOrder])
CompleteData = paste(as.character(CompleteData),'/38',sep="")
# Including a bottom line indicating how many valid subjects included for each contrast
ValidPart = rbind(ValidPart,CompleteData)
rownames(ValidPart)[length(rownames(ValidPart))] = "Valid subjects"
# Print table
kable(ValidPart,digits = 2, caption = "Index of the participant(s) rejected based on the following criteria:", col.names = rep("",13)) %>%
add_header_above(c("","Session 3" = 1, "GD-TD" = 1, "GC-TD" = 1, "GD-TC" = 1, "GD-TD" = 1, "GC-TD" = 1, "GD-TC" = 1, "GD-TD" = 1, "GC-TD" = 1, "GD-TC" = 1, "GD-TD" = 1, "GC-TD" = 1, "GD-TC" = 1)) %>%
add_header_above(c("","DPEQ" = 1, "BLOCGNG" = 3, "BLOCRVIP" = 3, "GNG_CR" = 3,  "RVIPN1" = 3)) %>%
add_header_above(c("","OVERALL" = 1, "BEHAVIOUR" = 6, "ERP" = 6)) %>%
kable_styling(bootstrap_options = c("striped", "hover"))
# Export the list of valid participants for each variable to be loaded in Analyses script
save(ValidPart,file=paste(RelativePath,"ValidPart.RData",sep="/"))
# Restricting to GD/TD !
NormalityData = RVIPData[RVIPData$Conditions=="GD-TD",c(3,6,17)]
# 1) Accurary (% HIT)
PlotViolinHline(NormalityData,NormalityData$Conditions,NormalityData$`%HIT`,"","RVIP HIT (in %)", 30)
# 2) Reaction Times (RT in ms)
PlotViolinHline(NormalityData,NormalityData$Conditions,NormalityData$meanRT_clean,"","RVIP RT (in ms)", c(350,750))
# Normality test
NormalityTBL = NormalityData %>%
normality()
kable(NormalityTBL,digits = 2, caption = "RVIP: Normality test (shapiro-Wilk)") %>%   kable_styling(bootstrap_options = c("striped", "hover"))
# One-sample t-tests
# 1) Accurary (% HIT)
RVIPHITT = t.test(NormalityData[,2], alternative = "greater", mu = 0.3)
TabOut = data.frame(t = as.numeric(RVIPHITT$statistic), Df = as.numeric(RVIPHITT$parameter), p = RVIPHITT$p.value)
kable(TabOut, caption = "RVIP HIT: One-sample t-test (>30%)") %>%
kable_styling(bootstrap_options = c("striped", "hover"))
# 2) Reaction Times (RT in ms)
RVIPRTT_low = t.test(NormalityData[,1], alternative = "greater", mu = 350)
TabOut = data.frame(t = as.numeric(RVIPRTT_low$statistic), Df = as.numeric(RVIPRTT_low$parameter), p = RVIPRTT_low$p.value)
kable(TabOut, caption = "RVIP RT: One-sample t-test (>350ms)") %>%
kable_styling(bootstrap_options = c("striped", "hover"))
RVIPRTT_high = t.test(NormalityData[,1], alternative = "less", mu = 750)
TabOut = data.frame(t = as.numeric(RVIPRTT_high$statistic), Df = as.numeric(RVIPRTT_high$parameter), p = RVIPRTT_high$p.value)
kable(TabOut, caption = "RVIP RT: One-sample t-test (<750ms)") %>%
kable_styling(bootstrap_options = c("striped", "hover"))
# Restricting to GD/TD !
NormalityData = GNGData[GNGData$Conditions=="GD-TD",c(3,9,20)]
# 1) False Alarms (% FA)
PlotViolinHline(NormalityData,NormalityData$Conditions,NormalityData$`%FA`,"","GNG FA (in %)", 10)
# 2) Reaction Times (RT in ms)
PlotViolinHline(NormalityData,NormalityData$Conditions,NormalityData$meanRT_clean,"","GNG RT (in ms)", c(170,570))
# Normality test
NormalityTBL = NormalityData %>%
normality()
kable(NormalityTBL,digits = 2, caption = "GNG: Normality test (shapiro-Wilk)") %>%   kable_styling(bootstrap_options = c("striped", "hover"))
# One-sample t-tests
# 1) False Alarms (% FA)
GNGFA = t.test(NormalityData[,2], alternative = "greater", mu = 10)
TabOut = data.frame(t = as.numeric(GNGFA$statistic), Df = as.numeric(GNGFA$parameter), p = GNGFA$p.value)
kable(TabOut, caption = "GNG FA: One-sample t-test (>10%)") %>%
kable_styling(bootstrap_options = c("striped", "hover"))
# 2) Reaction Times (RT in ms)
GNGRTT_low = t.test(NormalityData[,1], alternative = "greater", mu = 170)
TabOut = data.frame(t = as.numeric(GNGRTT_low$statistic), Df = as.numeric(GNGRTT_low$parameter), p = GNGRTT_low$p.value)
kable(TabOut, caption = "GNG RT: One-sample t-test (>170ms)") %>%
kable_styling(bootstrap_options = c("striped", "hover"))
GNGRTT_high = t.test(NormalityData[,1], alternative = "less", mu = 570)
TabOut = data.frame(t = as.numeric(GNGRTT_high$statistic), Df = as.numeric(GNGRTT_high$parameter), p = GNGRTT_high$p.value)
kable(TabOut, caption = "GNG HIT: One-sample t-test (<570ms)") %>%
kable_styling(bootstrap_options = c("striped", "hover"))
BCAFFData = data.frame (Data = c(Data$BCAFF_Withdrawal,Data$BCAFF_Energy,Data$BCAFF_Physical),
Label = rep(c("Withdrawal","Energy","Physical"),each = dim(Data)[1]))
BCAFFData = BCAFFData[!(is.na(BCAFFData[,1])),] # Remove NA lines
# Normality test
NormalityData = Data[,c(48,49,52,63)]
AvgData = data.frame (Conditions = NormalityData$Conditions,
BCaffMean = rowMeans(NormalityData[,-4], na.rm = T)) # Average over the 3 subscales
NormalityTBL = AvgData %>% normality()
kable(NormalityTBL,digits = 2, caption = "B-CAFF: Normality test (shapiro-Wilk)") %>%
kable_styling(bootstrap_options = c("striped", "hover"))
# Plot
PlotViolinHline(AvgData,"Average over 3 subscales",AvgData$BCaffMean,"","BCAFF Outcome Neutral Control", 3)
# One-sample t-tests
# 1) Withdrawal
BCAFFTest = t.test(AvgData$BCaffMean, alternative = "greater", mu = 3)
TabOut = data.frame(t = as.numeric(BCAFFTest$statistic), Df = as.numeric(BCAFFTest$parameter), p = BCAFFTest$p.value)
kable(TabOut, caption = "B-CAFF subscales average: One-sample t-test (mean score>3)") %>%
kable_styling(bootstrap_options = c("striped", "hover"))
load("F:/BIDS_Results/code/Diagnostics/Data/DiagnosticsDATA.RData")
TxtS
delete(Txts)
TxtS
delete(TxtS)
save.image("F:/BIDS_Results/code/Diagnostics/Data/DiagnosticsDATA.RData")
load("F:/BIDS_Results/code/Diagnostics/Data/DiagnosticsDATA.RData")
TxtS
delete(TxtS)
TxtS
remove(TxtS)
save.image("F:/BIDS_Results/code/Diagnostics/Data/DiagnosticsDATA.RData")
load("F:/BIDS_Results/code/Diagnostics/Data/DiagnosticsDATA.RData")
TxtS
