all(ParamCond[3],ParamCond[2]), # Contrast 2: "GD-TD" Vs "GD-TC"
all(ParamCond[2],ParamCond[1])) # Contrast 3: "GD-TC" Vs "GC-TD"
return(Param)
}
# Function for computing dependent-sample t-test or non-parametric equivalent
ComputeT = function(DF,VDStr,VIStr,Param,Comparisons,Alternate,Label=NULL,Paired=T){
# Preallocation
OUTTable = vector(mode = "list", length = length(Comparisons))
OUTStats = vector(mode = "list", length = length(Comparisons))
Pvals = rep(NA,length(Comparisons))
Formula = VDStr ~ VIStr
for (k in seq(1,length(Comparisons))) {
# Restrict Data to specific contrast
Idx = DF[,VIStr]==Comparisons[[k]][1] | DF[,VIStr]==Comparisons[[k]][2]
# Bayes Factor
# https://cran.r-project.org/web/packages/BayesFactor/vignettes/manual.html
DF1NARM = DF[DF[,VIStr]==Comparisons[[k]][1] & complete.cases(DF),]
DF2NARM = DF[DF[,VIStr]==Comparisons[[k]][2] & complete.cases(DF),]
# Only keeping subjects with data in both conditions (since paired t-test)
if (Paired == T) {
Idx1 = DF1NARM$Subjects %in% DF2NARM$Subjects
Idx2 = DF2NARM$Subjects %in% DF1NARM$Subjects
if (any(Idx1) | any(Idx2)) {
DF1NARM = DF1NARM[Idx1,]
DF2NARM = DF2NARM[Idx2,]
}
}
# Grouping both NA.rm df together
DFTEMP = rbind(DF1NARM,DF2NARM)
VDTEMP = DFTEMP[, VDStr]
VD1 = DF1NARM[,VDStr]
VD2 = DF2NARM[,VDStr]
Formula = as.formula(paste0(VDStr , " ~ ", VIStr))
# Confidence interval for one-sided test
# See: https://cran.r-project.org/web/packages/BayesFactor/vignettes/manual.html
# one-sided hypotheses that δ<0 versus the point null. We set nullInterval to c(-Inf,0):
NullInterv = c(-Inf,Inf)
if (Alternate[k] == "less") {NullInterv = c(-Inf,0) # case where x has a smaller mean than y
} else if (Alternate[k] == "greater") {NullInterv = c(0,Inf)} # case where x has a larger mean than y
# Compute, interpret and export BF
BayesF = ttestBF(DF1NARM[,VDStr],
DF2NARM[,VDStr],
paired = Paired, nullInterval = NullInterv) # one-sided t-test
BayesF = extractBF(BayesF)
BF01 = 1/BayesF$bf[1]
InterpretBF01 = interpret_bf(BF01, include_value = F)
InterpretBF01 = paste(InterpretBF01,"H0")
# First line is the probability that the mean difference be in the [-Inf 0] interval (H1)
# The second line is the probability that the mean difference is not in the interval (H0)
# Hence, BF01 = 0.1548 / 0.1846 = 0.84
# Here, the null hypothesis is 0.84 time(s) more probable than H1.
# https://cran.r-project.org/web/packages/bayestestR/vignettes/bayes_factors.html
# Computing statistics
if (Param[k]==T) {
#  paired sample t-test
Tbl = t.test(VD1,VD2,paired=Paired,alternative = Alternate[k]) # alternative = "greater" is the alternative that x has a larger mean than y
# Effect size (removing unused levels if factor type)
if (is.factor(DFTEMP[,VIStr])==T) {
ES_T = hedges_g(DFTEMP[,VDStr],droplevels(DFTEMP[,VIStr]),paired = Paired, alternative=Alternate[k], na.rm = T)
} else {
ES_T = hedges_g(DFTEMP[,VDStr],as.factor(DFTEMP[,VIStr]),paired = Paired, alternative=Alternate[k], na.rm = T)}
# OUTPUTS
# Statistics
OUTStats[[k]] = cbind(data.frame(t = as.numeric(Tbl$statistic), Df = as.numeric(Tbl$parameter), p_FDR = Tbl$p.value, HedgesG = ES_T$Hedges_g[1]), "2.5CI"=Tbl$conf.int[1], "97.5CI"=Tbl$conf.int[2], BayesFactor01 = BF01, InterpBF = as.character(InterpretBF01))
# Non-parametric statistics
} else if (Param[k]==F){
# Wilcoxon signed rank test
Tbl = wilcox.test(VD1,VD2,conf.int = T,paired = Paired,alternative  = Alternate[k])
# Effect size
if (Paired==T) {
ES_U = wilcoxonPairedR(x = VDTEMP,g = ifelse(DFTEMP[,VIStr]==Comparisons[[k]][1],0,1))}
else if (Paired==F) {
ES_U = wilcoxonR(x = VDTEMP,g = ifelse(DFTEMP[,VIStr]==Comparisons[[k]][1],0,1))}
# OUTPUTS
# Statistics
OUTStats[[k]] = cbind(data.frame(Stat = as.numeric(Tbl$statistic), p_FDR = Tbl$p.value, z = qnorm(Tbl$p.value/ 2), r = ES_U[1]), "2.5CI"=Tbl$conf.int[1], "97.5CI"=Tbl$conf.int[2], BayesFactor01 = BF01, InterpBF = as.character(InterpretBF01));
rownames(OUTStats[[k]]) = NULL
}
}
# Adjust p-values with FDR-correction
for (k in seq(1,length(Comparisons))) {
Pvals[k] = OUTStats[[k]]$p_FDR
}
Pvals = p.adjust(Pvals,"fdr")
# Replace p-values in list
for (k in seq(1,length(Comparisons))) {
OUTStats[[k]]$p_FDR = Pvals[k]
}
for (k in seq(1,length(Comparisons))) {
# Test name
if (Paired==T) {
if (Param[k]==T) {Test = "Paired samples t-test: "}
else {Test = "Wilcoxon signed rank test: "}
} else if (Paired==F) {
if (Param[k]==T) {Test = "Independent samples t-test: "}
else {Test = "Mann-Whitney U test: "}}
Caption = paste(Label,Test,Comparisons[[k]][1],ifelse(Alternate[k]=="less","<",">"),Comparisons[[k]][2])
# Build tables
kable(OUTStats[[k]], digits = 3, caption = Caption) %>%
kable_styling(bootstrap_options = c("striped", "hover")) %>%
print()
}
OUT = list(Table = OUTTable, Stats = OUTStats)
return(OUT)
}
# Function for computing permutation-based dependent-sample t-tests (for ERP data)
ComputeTPerm = function (DF,VDStr,VIStr,Comparisons,Alternate){
# Initialize output variable
OUTStats = data.frame(Contrast = rep(NA,length(Comparisons)), Perm_p = rep(NA,length(Comparisons)),
"Perm_2.5CI" = rep(NA,length(Comparisons)), "Perm_97.5CI" = rep(NA,length(Comparisons)),
Permutations = rep(5000,3))
for (k in seq(1,length(Comparisons))) {
# Contrast
OUTStats$Contrast[k] = paste(Comparisons[[k]][1],ifelse(Alternate[k]=="less","<",">"),Comparisons[[k]][2])
# Restrict Data to specific contrast
DF1NARM = DF[DF[,VIStr]==Comparisons[[k]][1] & complete.cases(DF),]
DF2NARM = DF[DF[,VIStr]==Comparisons[[k]][2] & complete.cases(DF),]
# Missing subjects data (need to remove extra data because of paired = T)
Idx1 = DF1NARM$Subjects %in% DF2NARM$Subjects
Idx2 = DF2NARM$Subjects %in% DF1NARM$Subjects
if (any(Idx1) | any(Idx2)) {
DF1NARM = DF1NARM[Idx1,]
DF2NARM = DF2NARM[Idx2,]
}
# Compute dependent-samples t-test with 5000 permutations
Tbl = perm.t.test(DF1NARM[,VDStr],DF2NARM[,VDStr], alternative = Alternate[k],paired = T, R = 5000)
# Saving data in output data frame
OUTStats$Perm_p[k] = Tbl$perm.p.value
OUTStats$Perm_2.5CI[k] = as.numeric(Tbl$perm.conf.int[1])
OUTStats$Perm_97.5CI[k] = as.numeric(Tbl$perm.conf.int[2])
OUTStats$Df = as.numeric(Tbl$parameter)
}
# Returning the dataframe
return(OUTStats)
}
# Loading the data into the workspace
load(paste(RelativePath,"Data/AnalysesDATA.RData",sep="/"))
# Compute descriptive statistics
Cols = c(4,6,10,13,17,18)
Li = seq(1,dim(QuestData)[1],3)
MeanDat = QuestData[Li,Cols] %>%  summarise_each(funs( mean( .,na.rm = TRUE)))
SDDat = QuestData[Li,Cols] %>%  summarise_each(funs( sd( .,na.rm = TRUE)))
MedianDat = QuestData[Li,Cols]  %>% summarise_each(funs( median( .,na.rm = TRUE)))
IQRDat = QuestData[Li,Cols] %>%  summarise_each(funs( IQR( .,na.rm = TRUE)))
# Build the dataframe
DescriptDF = data.frame(VD = colnames(MeanDat),
Mean = array(unlist(MeanDat)),
SD = array(unlist(SDDat)),
Median = array(unlist(MedianDat)),
IQR = array(unlist(IQRDat)))
# Print table
kable(DescriptDF,digits = 3, caption = "Descriptive statistics") %>%
kable_styling(bootstrap_options = c("striped", "hover"))
# Extract percentage of gender representation
PercData = QuestData[Li,] %>% dplyr::count(Gender) %>%
mutate(pct = prop.table(n))
GenderPlot = ggplot(PercData, aes(x = Gender, y = pct, label = scales::percent(pct,accuracy = 0.01))) +
geom_col(position = 'dodge') +
geom_text(position = position_dodge(width = .9), vjust = 1.5, size = 4) +
scale_y_continuous(labels = scales::percent) +
theme_minimal() +
theme(text = element_text(size = 14),
plot.title = element_text(hjust = 0.5),
legend.title = element_blank(),
plot.caption = element_text(hjust = 0)) +
labs(caption = "Note: Gender distribution", x = "Gender", y = "[%]")
# Extract percentage of gender representation
PercData = QuestData[Li,] %>% dplyr::count(DrugsHabits) %>%
mutate(pct = prop.table(n))
DrugsPlot = ggplot(PercData, aes(x = DrugsHabits, y = pct, label = scales::percent(pct,accuracy = 0.01))) +
geom_col(position = 'dodge') +
geom_text(position = position_dodge(width = .9), vjust = 1.5, size = 4) +
scale_y_continuous(labels = scales::percent) +
theme_minimal() +
theme(text = element_text(size = 14),
plot.title = element_text(hjust = 0.5),
legend.title = element_blank(),
plot.caption = element_text(hjust = 0)) +
labs(caption = "Note: Drugs habits distribution", x = "Drugs habits", y = "[%]")
# Extract percentage of gender representation
PercData = QuestData[Li,] %>% dplyr::count(WeeklySport) %>%
mutate(pct = prop.table(n))
SportPlot = ggplot(PercData, aes(x = WeeklySport, y = pct, label = scales::percent(pct,accuracy = 0.01))) +
geom_col(position = 'dodge') +
geom_text(position = position_dodge(width = .9), vjust = 1.5, size = 4) +
scale_y_continuous(labels = scales::percent) +
theme_minimal() +
theme(text = element_text(size = 14),
plot.title = element_text(hjust = 0.5),
legend.title = element_blank(),
plot.caption = element_text(hjust = 0)) +
labs(caption = "Note: Sport weekly distribution", x = "Sport weekly", y = "[%]")
# Plot alltogether
Graph = grid.arrange(GenderPlot, DrugsPlot, SportPlot, nrow = 1)
require(grid)
# PNG
grid.newpage()
# Open the device
png(paste(getwd(),"/FIGURES/FlowChart.png", sep=""), width = 200, height = 200, units='mm', res = 600)
# set some parameters to use repeatedly
midx = .25; rightx = .65; width = .4
gp = gpar(fill = "lightgrey")
# create boxes
(total = boxGrob("Initial sample\n N = 55", x=midx, y=.75, box_gp = gp, width = width))
(rando = boxGrob("Final sample\n N = 38 - 48", x=midx, y=.25, box_gp = gp, width = width))
# connect boxes like this
connectGrob(total, rando, "v")
# Removed participants
(inel = boxGrob("Dropouts:\n - Exclusion criteria (N = 1)\n - Side effects after first session (N = 2)\n - Unknown/personal reasons (N = 4)",
x=rightx, y=.58, box_gp = gp, width = .5, just = "left")) #, height = .05
connectGrob(total, inel, "-")
(inel = boxGrob("Exclusion based on outcome neutral controls\n (N = 0 - 10, depending on the contrast and DV)",
x=rightx, y=.42, box_gp = gp, width = .5, just = "left")) #, height = .05
connectGrob(total, inel, "-")
# close the device
Garbage = dev.off()
# removing Subjects column
AnalysisDataNoSubj = subset(AnalysisData, select=-Subjects)
# Computes stats
MeanDat = AnalysisDataNoSubj %>%  group_by(Conditions) %>% summarise_each(funs( mean( .,na.rm = TRUE)))
SDDat = AnalysisDataNoSubj %>%  group_by(Conditions) %>% summarise_each(funs( sd( .,na.rm = TRUE)))
MedianDat = AnalysisDataNoSubj %>%  group_by(Conditions) %>% summarise_each(funs( median( .,na.rm = TRUE)))
IQRDat = AnalysisDataNoSubj %>%  group_by(Conditions) %>% summarise_each(funs( IQR( .,na.rm = TRUE)))
# Computing N
CountDat = data.frame(Labels = levels(AnalysisDataNoSubj$Conditions))
CountDat = AnalysisDataNoSubj %>% group_by(Conditions) %>% summarise_all(funs(sum(!is.na(.))))
# Build the dataframe
Names = colnames(MeanDat)[2:dim(MeanDat)[2]]
DescriptDF = data.frame(VD = rep(Names,each=length(levels(AnalysisDataNoSubj$Conditions))),
VI = rep(levels(AnalysisDataNoSubj$Conditions),length(Names)),
Mean = array(unlist(MeanDat[,2:dim(MeanDat)[2]])),
SD = array(unlist(SDDat[,2:dim(MeanDat)[2]])),
Median = array(unlist(MedianDat[,2:dim(MeanDat)[2]])),
IQR = array(unlist(IQRDat[,2:dim(MeanDat)[2]])),
n = as.numeric(unlist(CountDat[,2:dim(MeanDat)[2]])))
# Print
kable(DescriptDF,digits = 3, caption = "Descriptive statistics (parametric and non-parametric)") %>%
kable_styling(bootstrap_options = c("striped", "hover"))
# Compute normality (Shapiro-Wilk test)
NormalityTBL = AnalysisDataNoSubj %>%
group_by(Conditions) %>%
normality()
# Adding skewness/kurtosis in the Normality table
DescriptStats = describeBy(AnalysisDataNoSubj, AnalysisDataNoSubj$Conditions, mat=TRUE, na.rm = TRUE, type=2)
# Type 2 estimate for skew/kurt = same as SPSS
DescriptStats = DescriptStats[-seq(1,3),] # Remove first three lines
NormalityTBL = cbind(NormalityTBL,DescriptStats$skew,DescriptStats$kurtosis)
colnames(NormalityTBL)[seq(6,7)] = c("Skewness","Kurtosis")
NormalityTBL$sample = DescriptStats$n # Replace the N-values from normality function since it counts also NA values
colnames(NormalityTBL)[1] = "Variable"
# Print
kable(NormalityTBL,digits = 3, caption = "Normality test (shapiro-Wilk) & Skewness/Kurtosis") %>%
kable_styling(bootstrap_options = c("striped", "hover"))
# VISUAL METHOD
# Density plots
# http://www.sthda.com/english/wiki/normality-test-in-r
GGD1 = ggdensity(AnalysisDataNoSubj$RVIPHIT, main = "Density plot of RVIPHIT",xlab = "RVIPHIT [%]")
GGD2 = ggdensity(AnalysisDataNoSubj$RVIPMeanRT, main = "Density plot of RVIPMeanRT",xlab = "RVIPMeanRT [ms]")
GGD3 = ggdensity(AnalysisDataNoSubj$GNGFA, main = "Density plot of GNGFA",xlab = "GNGFA [%]")
GGD4 = ggdensity(AnalysisDataNoSubj$GNGMeanRT, main = "Density plot of GNGMeanRT",xlab = "GNGMeanRT [ms]")
ggarrange(GGD1,GGD2,GGD3,GGD4,labels = c("A","B","C","D"),ncol = 2,nrow = 2)
# Q-Q plots
QQ1 = ggqqplot(AnalysisDataNoSubj$RVIPHIT, title = "Q-Q plot of RVIPHIT")
QQ2 = ggqqplot(AnalysisDataNoSubj$RVIPMeanRT, title = "Q-Q plot of RVIPMeanRT")
QQ3 = ggqqplot(AnalysisDataNoSubj$GNGFA, title = "Q-Q plot of GNGFA")
QQ4 = ggqqplot(AnalysisDataNoSubj$GNGMeanRT, title = "Q-Q plot of GNGMeanRT")
ggarrange(QQ1,QQ2,QQ3,QQ4,labels = c("A","B","C","D"),ncol = 2,nrow = 2)
# Parameters set to detect the presence of training effects across sessions
# Include a session column to the AnalysisData dataframe
AnalysisData$Sessions = CondData$SESSIONS
AnalysisDataNoSubj$Sessions = AnalysisData$Sessions
AnalysisDataNoSubjSess = subset(AnalysisDataNoSubj, select = -c(Conditions))
# Compute normality (Shapiro-Wilk test)
NormalityTBLSess = AnalysisDataNoSubjSess %>%
group_by(Sessions) %>%
normality()
# Adding skewness/kurtosis in the Normality table
DescriptStats = describeBy(AnalysisDataNoSubjSess, AnalysisDataNoSubjSess$Sessions, mat=TRUE, na.rm = TRUE, type=2)
# Type 2 estimate for skew/kurt = same as SPSS
DescriptStats = DescriptStats[-seq(13,15),] # Remove three last lines
NormalityTBLSess = cbind(NormalityTBLSess,DescriptStats$skew,DescriptStats$kurtosis)
colnames(NormalityTBLSess)[seq(6,7)] = c("Skewness","Kurtosis")
NormalityTBLSess$sample = DescriptStats$n # Replace the N-values from normality function since it counts also NA values
colnames(NormalityTBLSess)[1] = "Variable"
#### SAME but only for S1
# Compute normality (Shapiro-Wilk test)
AnalysisDataNoSubjSess1 = AnalysisDataNoSubj[AnalysisDataNoSubj$Sessions=="S1",]
NormalityTBLSess1 = AnalysisDataNoSubjSess1 %>%
group_by(Conditions) %>%
normality()
# Adding skewness/kurtosis in the Normality table
DescriptStats = describeBy(AnalysisDataNoSubjSess1, AnalysisDataNoSubjSess1$Conditions, mat=TRUE, na.rm = TRUE, type=2)
# Type 2 estimate for skew/kurt = same as SPSS
DescriptStats = DescriptStats[-c(1,2,3,16,17,18),] # Remove three last lines
NormalityTBLSess1 = cbind(NormalityTBLSess1,DescriptStats$skew,DescriptStats$kurtosis)
colnames(NormalityTBLSess1)[seq(6,7)] = c("Skewness","Kurtosis")
NormalityTBLSess1$sample = DescriptStats$n # Replace the N-values from normality function since it counts also NA values
colnames(NormalityTBLSess1)[1] = "Variable"
# Restricting AnalysisData dataframe to S1
AnalysisDataS1 = AnalysisData[AnalysisData$Sessions=="S1",]
# Check normality of distribution
Param = NormalityCheck(NormalityTBL,"RVIPHIT")
# Define alternate hypothesis
Alternate = c("less","less","less")
# Compute test statistics
Stats = ComputeT(AnalysisData,"RVIPHIT","Conditions",Param,Comparisons,Alternate)
# Plotting results
PlotViolin(AnalysisData,AnalysisData$Conditions,AnalysisData$RVIPHIT,"[%]","RVIP HIT",Stats,Comparisons)
# Check normality of distribution
Param = NormalityCheck(NormalityTBL,"RVIPMeanRT")
# Define alternate hypothesis
Alternate = c("greater","greater","greater")
# Compute test statistics
Stats = ComputeT(AnalysisData,"RVIPMeanRT","Conditions",Param,Comparisons,Alternate)
# Plotting results
PlotViolin(AnalysisData,AnalysisData$Conditions,AnalysisData$RVIPMeanRT,"[ms]","RVIP Reaction Times",Stats,Comparisons)
# Check normality of distribution
Param = NormalityCheck(NormalityTBL,"GNGFA")
# Define alternate hypothesis
Alternate = c("greater","greater","greater")
# Compute test statistics
Stats = ComputeT(AnalysisData,"GNGFA","Conditions",Param,Comparisons,Alternate)
# Plotting results
PlotViolin(AnalysisData,AnalysisData$Conditions,AnalysisData$GNGFA,"[%]","Go/NoGo FA",Stats,Comparisons)
# Check normality of distribution
Param = NormalityCheck(NormalityTBL,"GNGMeanRT")
# Define alternate hypothesis
Alternate = c("greater","greater","greater")
# Compute test statistics
Stats = ComputeT(AnalysisData,"GNGMeanRT","Conditions",Param,Comparisons,Alternate)
# Plotting results
PlotViolin(AnalysisData,AnalysisData$Conditions,AnalysisData$GNGMeanRT,"[ms]","Go/NoGo Reaction Times",Stats,Comparisons)
# Define alternate hypothesis
Alternate = c("less","less","less") # Same for all contrasts
# THERE MAY BE A MISTAKE AS THE N is not the same for GD-TC GNG (N2 & P3)
# between the ERP values (N = 44) and the source values (N = 43)
# The discrepancy relates to missing P900_S2 file !!
# The problem comes from errors in the allocation of files for analysis of source localization!!! e.g. P900_S3 is not in the right folder!!!
### ERP DATA
# Grouping the datasets
DescriptERP = rbind(cbind(RVIPDataP3[,c(2,5,7)]),GNGDataN2[,c(2,5,7)],GNGDataP3[,c(2,5,7)])
# Adding column specify the task and erp component
DescriptERP$Task = c(rep("RVIP",dim(RVIPDataP3)[1]),rep("GNG",dim(GNGDataN2)[1]+dim(GNGDataP3)[1]))
DescriptERP$ERPComp = c(rep("P3",dim(RVIPDataP3)[1]),rep("N2",dim(GNGDataN2)[1]),rep("P3",dim(GNGDataP3)[1]))
DescriptERP = DescriptERP[,c(4,5,1:3)]# Changing column order
# Turning table from wide to long
# Would benefit from a long format
DescriptERP = DescriptERP %>% gather(IndivGFPMeanOverPOI, IndivVoltageMeanOverPOIClust, -c(Task, ERPComp, Levels))
colnames(DescriptERP)[3:5] = c("Conditions","Variable","Value")
# Computes stats
MeanDat = DescriptERP %>%  group_by(Conditions,Task,ERPComp,Variable) %>% summarise_each(funs( mean( .,na.rm = TRUE)))
SDDat = DescriptERP %>%  group_by(Conditions,Task,ERPComp,Variable) %>% summarise_each(funs( sd( .,na.rm = TRUE)))
MedianDat = DescriptERP %>%  group_by(Conditions,Task,ERPComp,Variable) %>% summarise_each(funs( median( .,na.rm = TRUE)))
IQRDat = DescriptERP %>%  group_by(Conditions,Task,ERPComp,Variable) %>% summarise_each(funs( IQR( .,na.rm = TRUE)))
# Computing N
CountDat = DescriptERP %>% group_by(Conditions,Task,ERPComp,Variable) %>% summarise_all(funs(sum(!is.na(.))))
# # Build the dataframe
DescriptDF = MeanDat; colnames(DescriptDF)[dim(DescriptDF)[2]] = "Mean"
DescriptDF$SD = SDDat$Value
DescriptDF$Median = MedianDat$Value
DescriptDF$IQR = IQRDat$Value
DescriptDF$N = CountDat$Value
DescriptDF$Variable = as.factor(DescriptDF$Variable)
levels(DescriptDF$Variable) = c("GFP","Voltage amplitude")
### SOURCE LOCALISATION DATA
# Preparing the dataset
RVIPDataP3_ISGath = RVIPDataP3_IS[,-2] %>% gather(ROIs, CSD, -c(Subjects, Levels))
GNGDataN2_ISGath = GNGDataN2_IS[,-2] %>% gather(ROIs, CSD, -c(Subjects, Levels))
GNGDataP3_ISGath = GNGDataP3_IS[,-2] %>% gather(ROIs, CSD, -c(Subjects, Levels))
DescriptIS = rbind(RVIPDataP3_ISGath,GNGDataN2_ISGath,GNGDataP3_ISGath)
# Adding task and ERPComp columns
DescriptIS$Task = c(rep("RVIP",dim(RVIPDataP3_ISGath)[1]),rep("GNG",dim(GNGDataN2_ISGath)[1]+dim(GNGDataP3_ISGath)[1]))
DescriptIS$ERPComp = c(rep("P3",dim(RVIPDataP3_ISGath)[1]),rep("N2",dim(GNGDataN2_ISGath)[1]),rep("P3",dim(GNGDataP3_ISGath)[1]))
# Computes stats
MeanDat = DescriptIS[,-1] %>%  group_by(Levels,Task,ERPComp,ROIs) %>% summarise_each(funs( mean( .,na.rm = TRUE)))
SDDat = DescriptIS[,-1] %>%  group_by(Levels,Task,ERPComp,ROIs) %>% summarise_each(funs( sd( .,na.rm = TRUE)))
MedianDat = DescriptIS[,-1] %>%  group_by(Levels,Task,ERPComp,ROIs) %>% summarise_each(funs( median( .,na.rm = TRUE)))
IQRDat = DescriptIS[,-1] %>%  group_by(Levels,Task,ERPComp,ROIs) %>% summarise_each(funs( IQR( .,na.rm = TRUE)))
# Computing N
CountDat = DescriptIS[,-1] %>% group_by(Levels,Task,ERPComp,ROIs) %>% summarise_all(funs(sum(!is.na(.))))
# Build temporary dataframe to concatenate with DescriptDF
TEMP = cbind(MeanDat,SD = SDDat$CSD,Median = MedianDat$CSD,IQR = IQRDat$CSD,N = CountDat$CSD)
colnames(TEMP)[c(1,4,5)]=c("Conditions","Variable","Mean")
DescriptDF = rbind(DescriptDF,TEMP)
# Ordering according to ERPComp and Task
DescriptDF = DescriptDF[stri_order(paste(DescriptDF$ERPComp,DescriptDF$Task,sep="_"), numeric = TRUE),]
# Print
kable(DescriptDF,digits = 5, caption = "Descriptive statistics (parametric and non-parametric)") %>%
kable_styling(bootstrap_options = c("striped", "hover"))
# Empty vectors
NormalityTBL = NA; Skew = NA; Kurt = NA
# Compute normality tests (Shapiro-Wilk)
NormalityTBL = rbind(NormalityTBL,RVIPDataP3[,c(2,5,7)] %>% group_by(Levels) %>% normality())
NormalityTBL = rbind(NormalityTBL,GNGDataN2[,c(2,5,7)] %>% group_by(Levels) %>% normality())
NormalityTBL = rbind(NormalityTBL,GNGDataP3[,c(2,5,7)]  %>% group_by(Levels) %>% normality())
NormalityTBL = rbind(NormalityTBL,RVIPDataP3_IS[,c(3:6)]  %>% group_by(Levels) %>% normality())
NormalityTBL = rbind(NormalityTBL,GNGDataN2_IS[,c(3:6)]  %>% group_by(Levels) %>% normality())
NormalityTBL = rbind(NormalityTBL,GNGDataP3_IS[,c(3:5)]  %>% group_by(Levels) %>% normality())
# Adding skewness/kurtosis in the Normality table
## ERP Level
# RVIP P3
DescriptStats = describeBy(RVIPDataP3[,c(5,7)] , RVIPDataP3$Levels, mat=TRUE, na.rm = TRUE, type=2)
Skew = c(Skew,DescriptStats$skew)
Kurt = c(Kurt,DescriptStats$kurtosis)
# GNG N2
DescriptStats = describeBy(GNGDataN2[,c(5,7)], GNGDataN2$Levels, mat=TRUE, na.rm = TRUE, type=2)
Skew = c(Skew,DescriptStats$skew)
Kurt = c(Kurt,DescriptStats$kurtosis)
# GNG P3
DescriptStats = describeBy(GNGDataP3[,c(5,7)] , GNGDataP3$Levels, mat=TRUE, na.rm = TRUE, type=2)
Skew = c(Skew,DescriptStats$skew)
Kurt = c(Kurt,DescriptStats$kurtosis)
## SOURCE Level
# RVIP P3
DescriptStats = describeBy(RVIPDataP3_IS[,c(4:6)] , RVIPDataP3_IS$Levels, mat=TRUE, na.rm = TRUE, type=2)
Skew = c(Skew,DescriptStats$skew)
Kurt = c(Kurt,DescriptStats$kurtosis)
# GNG N2
DescriptStats = describeBy(GNGDataN2_IS[,c(4:6)], GNGDataN2_IS$Levels, mat=TRUE, na.rm = TRUE, type=2)
Skew = c(Skew,DescriptStats$skew)
Kurt = c(Kurt,DescriptStats$kurtosis)
# GNG P3
DescriptStats = describeBy(GNGDataP3_IS[,c(4,5)] , GNGDataP3_IS$Levels, mat=TRUE, na.rm = TRUE, type=2)
Skew = c(Skew,DescriptStats$skew)
Kurt = c(Kurt,DescriptStats$kurtosis)
# Adding Skewness/Kurtosis columns
NormalityTBL$Skewness = Skew
NormalityTBL$Kurtosis = Kurt
NormalityTBL = NormalityTBL[-1,] # removing first line (all NAs)
# Adding column specify the task and erp component
NormalityTBL$Task = c(rep("RVIP",6),rep("GNG",12),rep("RVIP",9),rep("GNG",9+6))
NormalityTBL$ERPComp = c(rep("P3",6),rep("N2",6),rep("P3",6),rep("P3",9),rep("N2",9),rep("P3",6))
NormalityTBL = NormalityTBL[,c(2,8,9,1,3:7)]# Changing column order
# Renaming factor levels
NormalityTBL$variable = as.factor(NormalityTBL$variable)
Idx = c(which(levels(NormalityTBL$variable)=="IndivGFPMeanOverPOI"),which(levels(NormalityTBL$variable)=="IndivVoltageMeanOverPOIClust"))
levels(NormalityTBL$variable)[Idx] = c("GFP","Voltage amplitude")
colnames(NormalityTBL)[c(1,4)] = c("Conditions","Variable")
# Ordering according to ERPComp and Task
# NormalityTBL = NormalityTBL[stri_order(paste(NormalityTBL$Task,NormalityTBL$ERPComp,sep="_"), numeric = TRUE),]
# Table
kable(NormalityTBL,digits = 5, caption = "Normality test (shapiro-Wilk) & Skewness/Kurtosis") %>%
kable_styling(bootstrap_options = c("striped", "hover"))
# Check normality of distribution
Param = NormalityCheck(NormalityTBL,c("Voltage amplitude","RVIP","P3"))
# Run statistics
StatsVolt = ComputeT(RVIPDataP3,"IndivVoltageMeanOverPOIClust","Levels",Param,Comparisons,Alternate)
# Plots and tables
PlotViolin(RVIPDataP3,RVIPDataP3$Levels,RVIPDataP3$IndivVoltageMeanOverPOIClust,"[\U003BCV]","RVIP P3 mean voltage amplitude",StatsVolt,Comparisons)
# Run statistics
StatsGFP = ComputeTPerm(RVIPDataP3,"IndivGFPMeanOverPOI","Levels",Comparisons,Alternate)
# Print table
kable(StatsGFP,digits = 3, caption = "Permutation-based statistics for RVIP P3 GFP") %>%
kable_styling(bootstrap_options = c("striped", "hover"))
# Plots
PlotViolinPerm(RVIPDataP3,RVIPDataP3$Levels,RVIPDataP3$IndivGFPMeanOverPOI,"GFP","RVIPData P3 GFP",StatsGFP,Comparisons)
Plots = list()
for (k in 1:length(ROI_RVIPP3)) {
# Check normality of distribution
Param = NormalityCheck(NormalityTBL,c(ROI_RVIPP3[k],"RVIP","P3"))
# Run statistics
Stats[[k]] = ComputeT(RVIPDataP3_IS,ROI_RVIPP3[k],"Levels",Param,Comparisons,Alternate,Label = paste(toupper(ROI_RVIPP3[k]),":"))
# Plots and tables
Idx = which(colnames(RVIPDataP3_IS)==ROI_RVIPP3[k])
PlotViolin(RVIPDataP3_IS,RVIPDataP3_IS$Levels,RVIPDataP3_IS[,Idx],"Current source density (CSD)",ROI_RVIPP3[k],Stats[[k]],Comparisons)
}
# Check normality of distribution
Param = NormalityCheck(NormalityTBL,c("Voltage amplitude","GNG","N2"))
# Run statistics
StatsVolt = ComputeT(GNGDataN2,"IndivVoltageMeanOverPOIClust","Levels",Param,Comparisons,Alternate)
# Plots and tables
PlotViolin(GNGDataN2,GNGDataN2$Levels,GNGDataN2$IndivVoltageMeanOverPOIClust,"[\U003BCV]","Go/NoGo N2 mean voltage amplitude",StatsVolt,Comparisons)
# Run statistics
StatsGFP = ComputeTPerm(GNGDataN2,"IndivGFPMeanOverPOI","Levels",Comparisons,Alternate)
# Print table
kable(StatsGFP,digits = 3, caption = "Permutation-based statistics for GNG N2 GFP") %>%
kable_styling(bootstrap_options = c("striped", "hover"))
# Plots
PlotViolinPerm(GNGDataN2,GNGDataN2$Levels,GNGDataN2$IndivGFPMeanOverPOI,"GFP","Go/NoGo N2 GFP",StatsGFP,Comparisons)
Plots = list()
for (k in 1:length(ROI_GNGN2)) {
# Check normality of distribution
Param = NormalityCheck(NormalityTBL,c(ROI_GNGN2[k],"GNG","N2"))
# Run statistics
Stats[[k]] = ComputeT(GNGDataN2_IS,ROI_GNGN2[k],"Levels",Param,Comparisons,Alternate,Label = paste(toupper(ROI_GNGN2[k]),":"))
# Plots and tables
Idx = which(colnames(GNGDataN2_IS)==ROI_GNGN2[k])
PlotViolin(GNGDataN2_IS,GNGDataN2_IS$Levels,GNGDataN2_IS[,Idx],"Current source density (CSD)",ROI_GNGN2[k],Stats[[k]],Comparisons)
}
# Check normality of distribution
Param = NormalityCheck(NormalityTBL,c("Voltage amplitude","GNG","P3"))
# Run statistics
StatsVolt = ComputeT(GNGDataP3,"IndivVoltageMeanOverPOIClust","Levels",Param,Comparisons,Alternate)
# Plots and tables
PlotViolin(GNGDataP3,GNGDataP3$Levels,GNGDataP3$IndivVoltageMeanOverPOIClust,"[\U003BCV]","Go/NoGo P3 mean voltage amplitude",StatsVolt,Comparisons)
# Run statistics
StatsGFP = ComputeTPerm(GNGDataP3,"IndivGFPMeanOverPOI","Levels",Comparisons,Alternate)
# Print table
kable(StatsGFP,digits = 3, caption = "Permutation-based statistics for GNG P3 GFP") %>%
kable_styling(bootstrap_options = c("striped", "hover"))
# Plots
PlotViolinPerm(GNGDataP3,GNGDataP3$Levels,GNGDataP3$IndivGFPMeanOverPOI,"GFP","Go/NoGo P3 GFP",StatsGFP,Comparisons)
Plots = list()
for (k in 1:length(ROI_GNGP3)) {
# Check normality of distribution
Param = NormalityCheck(NormalityTBL,c(ROI_GNGP3[k],"GNG","P3"))
# Run statistics
Stats[[k]] = ComputeT(GNGDataP3_IS,ROI_GNGP3[k],"Levels",Param,Comparisons,Alternate,Label = paste(toupper(ROI_GNGP3[k]),":"))
# Plots and tables
Idx = which(colnames(GNGDataP3_IS)==ROI_GNGP3[k])
PlotViolin(GNGDataP3_IS,GNGDataP3_IS$Levels,GNGDataP3_IS[,Idx],"Current source density (CSD)",ROI_GNGP3[k],Stats[[k]],Comparisons)
}
