library(reshape2)
library(plyr)
library(lme4)
library(lmerTest)
library(ggplot2)
library(ggpubr)
library(emmeans)

# Read csv as data frame into environment - Note: change FILENAME
Data_long <- read.csv("C:/Users/jamessteele/Dropbox/Research Toolbox/Current R Data Files/All Data - long form.csv")

Data_long$Load <- as.character(Data_long$Load)
Data_long$Load <- revalue(Data_long$Load, c("1"="30% MVC", "2"="70% MVC"))
Data_long$Load <- as.factor(Data_long$Load)
Data_long$Condition <- as.character(Data_long$Condition)
Data_long$Condition <- revalue(Data_long$Condition, c("1"="Placebo", "2"="Caffeine"))
Data_long$Condition <- as.factor(Data_long$Condition)
Data_long$Participant <- as.factor(Data_long$Participant)
Data_long$effort <- as.integer(Data_long$effort)
Data_long$disc <- as.integer(Data_long$disc)
Data_long$Load.Condition <- paste(Data_long$Load, Data_long$Condition)
str(Data_long)

Data30 <- subset(Data_long, Load =="30% MVC")
Data70 <- subset(Data_long, Load =="70% MVC")

#Peak Torque
# Run lmm with random intercepts - lmer(dependent variable ~ Condition + Load + Condition:Load + (1|Participant), data = Data_long)
LMMPF <- lmer(PF ~ Condition + (1|Participant), data = Data_long, REML = TRUE)

# Calculate estimated marginal means - emmeans(LMMPF, "factor")
EMMEANSPF <- emmeans(LMMPF, "Condition")
EMMEANSPFdf <- as.data.frame(EMMEANSPF)
CONTRASTSPF <- emmeans(LMMPF, "Condition", contr = "pairwise")
CONTRASTSPFdf <- as.data.frame(CONTRASTSPF)

# Plot elements
TITLES <- element_text(size = 16)
AXIS <- element_text(size = 12)
BACKGROUND <- element_rect(fill = "White", colour = "Black")
THEME <- theme(axis.text = AXIS, axis.title = TITLES, panel.background = BACKGROUND,  
      legend.position = "top", legend.title = element_blank())

# Paired data plot
p <- plot(EMMEANSPF) + THEME + 
  labs(x = "Peak Torque (Nm)", y = "Condition") +
  geom_line(aes(x = PF, y = Condition, group = Participant), size = 0.5, colour = "indianred1", 
              data = Data30) +
  geom_line(aes(x = PF, y = Condition, group = Participant), size = 0.5, colour = "skyblue", 
            data = Data70) +
   coord_flip() +
  scale_x_continuous(limit = c(40, 200), breaks = c(40,60,80,100,120,140,160,180,200))

pairedplot <- p + geom_point(size = 1) + geom_path(aes(x = emmean, y = Condition, group = 1), size = 1, colour = "grey20", data = EMMEANSPFdf)

#Paired contrasts plot
#Paired contrasts with 90%CIs for equivalence
CONTRASTSPFequiv <- confint(CONTRASTSPF, level = .90)
CONTRASTSPFequivdf <- as.data.frame(CONTRASTSPFequiv)

p <- plot(CONTRASTSPF$contrasts) + THEME + 
  labs(x = "Difference in Peak Torque (Nm)", y = "Condition") +
  coord_flip() +
  scale_x_continuous(limit = c(-10, 25), breaks = c(-10,-5,0,5,10,15,20,25)) +
  geom_errorbarh(height=.05, size = 1, aes(xmin=CONTRASTSPFequiv$contrasts$lower.CL, xmax=CONTRASTSPFequiv$contrasts$upper.CL), colour="grey20")

contrastplot <- p + geom_vline(xintercept = 0, linetype = "dashed", color = "Black", size = 1)

# Add bands for SESOI - NOTE: 6.98% is halfwidth of minimal detectable change from reliability data in our lab for the IKD. 
# This determined the %contrast difference relative to the intercept from LMM 
# Then what 6.98% is in the raw units and add's these lines to the graph
CONTRASTSPFdf <- CONTRASTSPFdf[-c(1),] 
LMMdf <- as.data.frame(coef(summary(LMMPF))[ , "Estimate"])
CONTRASTSPFdf$percent <- (CONTRASTSPFdf$contrasts.estimate/LMMdf$`coef(summary(LMMPF))[, "Estimate"]`[1])*100
CONTRASTSPFdf$onepercent <- (CONTRASTSPFdf$contrasts.estimate/CONTRASTSPFdf$percent)
contrastplot <- contrastplot + geom_vline(xintercept = (CONTRASTSPFdf$onepercent * 6.98), linetype = "dashed", color = "Red", size = 1)
contrastplot <- contrastplot + geom_vline(xintercept = (CONTRASTSPFdf$onepercent * -6.98), linetype = "dashed", color = "Red", size = 1)

#Combine into one plot
ggarrange(pairedplot,contrastplot, ncol=2, nrow=1, labels = c("(A)","(B)"), 
          font.label = list(size = 10),common.legend = TRUE, legend = "bottom") 

#Time to Task Failure
# Run lmm with random intercepts - lmer(dependent variable ~ Condition + Load + Condition:Load + (1|Participant), data = Data_long)
LMMTF <- lmer(TF ~ Condition + Load + Condition:Load + (1|Participant), data = Data_long, REML = TRUE)


# Calculate estimated marginal means - emmeans(LMMTF, "factor")
EMMEANSTF <- emmeans(LMMTF, "Condition", by = "Load")
EMMEANSTFdf <- as.data.frame(EMMEANSTF)
CONTRASTSTF <- emmeans(LMMTF, "Condition", by = "Load", contr = "pairwise")
CONTRASTSTFdf <- as.data.frame(CONTRASTSTF)

# Paired data plot
#Paired contrasts with 90%CIs for equivalence
CONTRASTSTFequiv <- confint(CONTRASTSTF, level = .90)
CONTRASTSTFequivdf <- as.data.frame(CONTRASTSTFequiv)

p <- plot(EMMEANSTF) + THEME + 
  labs(x = "Time to Task Failure (seconds)", y = "Condition") +
  geom_line(aes(x = TF, y = Condition, group = Participant), size = 0.5, colour = "indianred1", 
            data = Data30) +
  geom_line(aes(x = TF, y = Condition, group = Participant), size = 0.5, colour = "skyblue", 
            data = Data70) +
  coord_flip() +
  scale_x_continuous(limit = c(0, 450), breaks = c(0,50,100,150,200,250,300,350,400,450))

pairedplot <- p + geom_point(size = 1) + geom_path(aes(x = emmean, y = Condition, group = 1), size = 1, colour = "grey20", data = EMMEANSTFdf)

#Paired contrasts plot
p <- plot(CONTRASTSTF$contrasts) + THEME + 
  labs(x = "Difference in Time to Task Failure (seconds)", y = "Condition") +
  coord_flip() +
  scale_x_continuous(limit = c(-40, 70), breaks = c(-40,-30,-20,-10,0,10,20,30,40,50,60,70)) +
  geom_errorbarh(height=.05, size = 1, aes(xmin=CONTRASTSTFequiv$contrasts$lower.CL, xmax=CONTRASTSTFequiv$contrasts$upper.CL), colour="grey20")

contrastplot <- p + geom_vline(xintercept = 0, linetype = "dashed", color = "Black", size = 1)

# Add bands for SESOI - NOTE: 6.98% is halfwidth of minimal detectable change from reliability data in our lab for the IKD. 
# This determined the %contrast difference relative to the intercept from LMM 
# Then what 6.98% is in the raw units and add's these lines to the graph
CONTRASTSTFdf <- CONTRASTSTFdf[-c(2,3),] 
LMMdf <- as.data.frame(coef(summary(LMMTF))[ , "Estimate"])
CONTRASTSTFdf$percent <- (CONTRASTSTFdf$contrasts.estimate/LMMdf$`coef(summary(LMMTF))[, "Estimate"]`[1])*100
CONTRASTSTFdf$onepercent <- (CONTRASTSTFdf$contrasts.estimate/CONTRASTSTFdf$percent)
contrastplot <- contrastplot + geom_vline(xintercept = (CONTRASTSTFdf$onepercent * 6.98), linetype = "dashed", color = "Red", size = 1)
contrastplot <- contrastplot + geom_vline(xintercept = (CONTRASTSTFdf$onepercent * -6.98), linetype = "dashed", color = "Red", size = 1)

#Combine into one plot
ggarrange(pairedplot,contrastplot, ncol=2, nrow=1, labels = c("(A)","(B)"), 
          font.label = list(size = 10),common.legend = TRUE, legend = "bottom") 

# Effort
# Run lmm with random intercepts - lmer(dependent variable ~ Condition + Load + Condition:Load + (1|Participant), data = Data_long)
LMMeffort <- lmer(effort ~ Condition + Load + Condition:Load + (1|Participant), data = Data_long, REML = TRUE)

# Calculate estimated marginal means - emmeans(LMMeffort, "factor")
EMMEANSeffort <- emmeans(LMMeffort, "Condition", by = "Load")
EMMEANSeffortdf <- as.data.frame(EMMEANSeffort)
CONTRASTSeffort <- emmeans(LMMeffort, "Condition", by = "Load", contr = "pairwise")
CONTRASTSeffortdf <- as.data.frame(CONTRASTSeffort)

# Paired data plot
#Paired contrasts with 90%CIs for equivalence
CONTRASTSeffortequiv <- confint(CONTRASTSeffort, level = .90)
CONTRASTSeffortequivdf <- as.data.frame(CONTRASTSeffortequiv)

p <- plot(EMMEANSeffort) + THEME + 
  labs(x = "Rating of Perceived Effort", y = "Condition") +
  geom_line(aes(x = effort, y = Condition, group = Participant), size = 0.5, colour = "indianred1", 
            data = Data30) +
  geom_line(aes(x = effort, y = Condition, group = Participant), size = 0.5, colour = "skyblue", 
            data = Data70) +
  coord_flip() +
  scale_x_continuous(limit = c(-0.5, 10.5), breaks = c(0,1,2,3,4,5,6,7,8,9,10))

pairedplot <- p + geom_point(size = 1) + geom_path(aes(x = emmean, y = Condition, group = 1), size = 1, colour = "grey20", data = EMMEANSeffortdf)

#Paired contrasts plot
p <- plot(CONTRASTSeffort$contrasts) + THEME + 
  labs(x = "Difference in Rating of Perceived Effort", y = "Condition") +
  coord_flip() +
  scale_x_continuous(limit = c(-2,2), breaks = c(-2,-1,0,1,2)) +
  geom_errorbarh(height=.05, size = 1, aes(xmin=CONTRASTSeffortequiv$contrasts$lower.CL, xmax=CONTRASTSeffortequiv$contrasts$upper.CL), colour="grey20")

contrastplot <- p + geom_vline(xintercept = 0, linetype = "dashed", color = "Black", size = 1)

# Add bands for SESOI - NOTE: Based on conservative of 0.85 TEM from http://trainology.org/PDF/v6-1%2001%20p1-8%20Steele%20st%20al%200419.pdf
# Halfwidth of minimal detecable change
contrastplot <- contrastplot + geom_vline(xintercept = ((0.85 * 1.96 * sqrt(2))/2), linetype = "dashed", color = "Red", size = 1)
contrastplot <- contrastplot + geom_vline(xintercept = ((-0.85 * 1.96 * sqrt(2))/2), linetype = "dashed", color = "Red", size = 1)

#Combine into one plot
ggarrange(pairedplot,contrastplot, ncol=2, nrow=1, labels = c("(A)","(B)"), 
          font.label = list(size = 10),common.legend = TRUE, legend = "bottom")

# Discomfort
# Run lmm with random intercepts - lmer(dependent variable ~ Condition + Load + Condition:Load + (1|Participant), data = Data_long)
LMMdisc <- lmer(disc ~ Condition + Load + Condition:Load + (1|Participant), data = Data_long, REML = TRUE)

# Calculate estimated marginal means - emmeans(LMMdisc, "factor")
EMMEANSdisc <- emmeans(LMMdisc, "Condition", by = "Load")
EMMEANSdiscdf <- as.data.frame(EMMEANSdisc)
CONTRASTSdisc <- emmeans(LMMdisc, "Condition", by = "Load", contr = "pairwise")
CONTRASTSdiscdf <- as.data.frame(CONTRASTSdisc)

# Paired data plot
#Paired contrasts with 90%CIs for equivalence
CONTRASTSdiscequiv <- confint(CONTRASTSdisc, level = .90)
CONTRASTSdiscequivdf <- as.data.frame(CONTRASTSdiscequiv)

p <- plot(EMMEANSdisc) + THEME + 
  labs(x = "Rating of Perceived Discomfort", y = "Condition") +
  geom_line(aes(x = disc, y = Condition, group = Participant), size = 0.5, colour = "indianred1", 
            data = Data30) +
  geom_line(aes(x = disc, y = Condition, group = Participant), size = 0.5, colour = "skyblue", 
            data = Data70) +
  coord_flip() +
  scale_x_continuous(limit = c(-0.5, 10.5), breaks = c(0,1,2,3,4,5,6,7,8,9,10))

pairedplot <- p + geom_point(size = 1) + geom_path(aes(x = emmean, y = Condition, group = 1), size = 1, colour = "grey20", data = EMMEANSdiscdf)

#Paired contrasts plot
p <- plot(CONTRASTSdisc$contrasts) + THEME + 
  labs(x = "Difference in Rating of Perceived Discomfort", y = "Condition") +
  coord_flip() +
  scale_x_continuous(limit = c(-2,2), breaks = c(-2,-1,0,1,2)) +
  geom_errorbarh(height=.05, size = 1, aes(xmin=CONTRASTSdiscequiv$contrasts$lower.CL, xmax=CONTRASTSdiscequiv$contrasts$upper.CL), colour="grey20")

contrastplot <- p + geom_vline(xintercept = 0, linetype = "dashed", color = "Black", size = 1)

# Add bands for SESOI - NOTE: Based on conservative of 1.01 TEM from http://trainology.org/PDF/v6-1%2001%20p1-8%20Steele%20st%20al%200419.pdf
# Halfwidth of minimal detecable change
contrastplot <- contrastplot + geom_vline(xintercept = ((1.01 * 1.96 * sqrt(2))/2), linetype = "dashed", color = "Red", size = 1)
contrastplot <- contrastplot + geom_vline(xintercept = ((-1.01 * 1.96 * sqrt(2))/2), linetype = "dashed", color = "Red", size = 1)

#Combine into one plot
ggarrange(pairedplot,contrastplot, ncol=2, nrow=1, labels = c("(A)","(B)"), 
          font.label = list(size = 10),common.legend = TRUE, legend = "bottom")

# Test Results and emmeans
print(LMMPF)
summary(LMMPF)
anova(LMMPF)
emmeans(LMMPF, "Condition", contr = "pairwise")
test(CONTRASTSPF)
deltaPF <- (CONTRASTSPFdf$onepercent * 6.98)
test(CONTRASTSPF, delta = deltaPF[1], side = "equivalence")
eff_size(CONTRASTSPF$emmeans, sigma = sigma(LMMPF), edf = 71)

print(LMMTF)
summary(LMMTF)
anova(LMMTF)
emmeans(LMMTF, "Condition", by = "Load", contr = "pairwise")
test(CONTRASTSTF)
deltaTF <- (CONTRASTSTFdf$onepercent * 6.98)
test(CONTRASTSTF, delta = deltaTF[1], side = "equivalence")
eff_size(CONTRASTSTF$emmeans, sigma = sigma(LMMTF), edf = 71)

print(LMMeffort)
summary(LMMeffort)
anova(LMMeffort)
emmeans(LMMeffort, "Condition", by = "Load", contr = "pairwise")
test(CONTRASTSeffort)
deltaeffort <- ((0.85 * 1.96 * sqrt(2))/2)
test(CONTRASTSeffort, delta = deltaeffort[1], side = "equivalence")
eff_size(CONTRASTSeffort$emmeans, sigma = sigma(LMMeffort), edf = 71)

print(LMMdisc)
summary(LMMdisc)
anova(LMMdisc)
emmeans(LMMdisc, "Condition", by = "Load", contr = "pairwise")
test(CONTRASTSdisc)
deltadisc <- ((1.01 * 1.96 * sqrt(2))/2)
test(CONTRASTSdisc, delta = deltadisc[1], side = "equivalence")
eff_size(CONTRASTSdisc$emmeans, sigma = sigma(LMMdisc), edf = 71)

# Q-Q plots for models
qqPF <- ggplot(LMMPF, aes(sample=PF))+stat_qq() + THEME + stat_qq_line(linetype = "dashed", color = "Red", size = 1)
qqTF <- ggplot(LMMTF, aes(sample=TF))+stat_qq() + THEME + stat_qq_line(linetype = "dashed", color = "Red", size = 1)
qqeffort <- ggplot(LMMeffort, aes(sample=effort))+stat_qq() + THEME + stat_qq_line(linetype = "dashed", color = "Red", size = 1)
qqdisc <- ggplot(LMMdisc, aes(sample=disc))+stat_qq() + THEME + stat_qq_line(linetype = "dashed", color = "Red", size = 1)

#Combine into one plot
ggarrange(qqPF,qqTF,qqeffort,qqdisc, ncol=2, nrow=2, labels = c("(A)","(B)","(C)","(D)"), font.label = list(size = 10))