JA_TMS_2022/Pilot_analysis/JA_pilot_analysis.R

rm(list = ls())


library(openxlsx)
library(dplyr) 
library(ggplot2)
library(tidyverse)
library(ez)
library(plotly)
library(ggsignif)
library(cowplot)
library(multipanelfigure)
library(ggh4x)
library(svglite)
library(BayesFactor)


setwd("/Users/b1037210/Desktop/Progetti/JA_Stage1/")

df = read.xlsx("final_table_for_R_10-2500_def_o1_rep.xlsx")


df$TMS_timings_rep = as.factor(df$TMS_timings_rep)
levels(df$TMS_timings_rep) = c("Late TMS","Early TMS")
df$TMS_timings_rep <- factor(df$TMS_timings_rep, levels = c("Early TMS","Late TMS"))

df$subject = as.factor(df$subject)
df$CF_colors_order = as.factor(df$CF_colors_order)
levels(df$CF_colors_order) = c("Lift", "Catch", "Press")

df <- subset(df, CF_colors_order == "Press"| CF_colors_order == "Lift" ) 

df$CF_colors_order = droplevels(df$CF_colors_order)

groups <- group_by(df, CF_colors_order, subject, TMS_timings_rep) 

# summarize by using the medians of each subject --> have a table with a "mean" column for each cell for each participant
data_avgs <- summarise(groups,
                       mean = median(ECD_MEP_amplitudes, na.rm=TRUE))

# take a peek
data_avgs = as.data.frame(data_avgs)

df_final = data_avgs

attach(df_final)


## descriptive statistics
s_alt = ezStats(data=df_final,
                dv = .(mean),
                wid = .(subject),
                within = .(CF_colors_order, TMS_timings_rep))

## plot now
pd = position_dodge(width=0.0)

# create basic plot with x axis as Cngruendy and y as mean (which is median actually)
g = ggplot(df_final, aes(x = CF_colors_order, y = mean))

# create Facets, the level of one variable above, the name of the levels of the other on the new line (\n)
g = g+facet_grid(col = vars(interaction(TMS_timings_rep, sep = "\n",lex.order = TRUE)))

# create bar plot 
g = g+ geom_bar(aes(fill = CF_colors_order), stat = "summary", fun = "mean")


g = g+ geom_line(aes(group = subject)
                 , stat = "summary"
                 , fun = "mean" 
                 ,size= .3
                 , position = pd)


#new filling colors for the bars --> values = c("red","orange","blue","green"))
g = g+scale_fill_manual(values = c("red","blue","tomato","dodgerblue"))

g = g+theme_classic()

g = g+theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) #top, right, bottom, left

g = g+theme(strip.text.x = element_text(size=20, angle=0, vjust = 0),
            strip.background = element_rect(colour="white", fill="white"))

g = g + theme(
  axis.title.x = element_text( size = 20, face = "bold", margin = margin(t = 10, r = 0, b = 0, l = 0)),
  axis.title.y = element_text(size = 20, face = "bold", margin = margin(t = 0, r = 10, b = 0, l = 0)),
  
  axis.text.x = element_text(size = 15),# family = 'Helvetica',
  axis.text.y = element_text(size = 15))  

g = g+xlab("Cue instruction") + ylab("peak-to-peak amplitude (µV)") #ylab("TMS-evoked pressure (a.u.)") #

g = g+theme(
  legend.position = "none")


file_name = "ECU" #"Pressure (z)"#"ECU(z) - FDS(z)"#"Flexor Digitorum Superficialis (Z-scores)\n"
g = g+ggtitle(file_name)
g = g+theme(plot.title = element_text(hjust = 0.5, size = 35, face = "bold"))


g


## stats
late_lift = subset(data_avgs, TMS_timings_rep == "Late TMS" & CF_colors_order == "Lift")
late_press = subset(data_avgs, TMS_timings_rep == "Late TMS" & CF_colors_order == "Press")

early_lift = subset(data_avgs, TMS_timings_rep == "Early TMS" & CF_colors_order == "Lift")
early_press = subset(data_avgs, TMS_timings_rep == "Early TMS" & CF_colors_order == "Press")

shapiro.test(late_lift$mean)
shapiro.test(late_press$mean)

shapiro.test(early_lift$mean)
shapiro.test(early_press$mean)

t.test(early_lift$mean, early_press$mean, paired = TRUE, alternative = "greater")
t.test(late_lift$mean, late_press$mean, paired = TRUE, alternative = "greater")



## save
file_title = paste0(file_name, ".jpg");
ggsave(file_title,g,width = 10, height = 7, device = "jpg")