###Packages

library(lme4)
## Warning: package 'lme4' was built under R version 4.1.1
## Loading required package: Matrix
library(readxl)
## Warning: package 'readxl' was built under R version 4.1.1
library(haven)
## Warning: package 'haven' was built under R version 4.1.1
library(effects)
## Warning: package 'effects' was built under R version 4.1.1
## Loading required package: carData
## Warning: package 'carData' was built under R version 4.1.1
## lattice theme set by effectsTheme()
## See ?effectsTheme for details.
library(lattice)
## Warning: package 'lattice' was built under R version 4.1.1
library(car)
## Warning: package 'car' was built under R version 4.1.1
library(ggplot2) 
## Warning: package 'ggplot2' was built under R version 4.1.1
library(knitr)
## Warning: package 'knitr' was built under R version 4.1.1
library(reshape2)
## Warning: package 'reshape2' was built under R version 4.1.1
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.1.1
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:car':
## 
##     recode
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(forcats)
## Warning: package 'forcats' was built under R version 4.1.1
library(DHARMa)
## Warning: package 'DHARMa' was built under R version 4.1.1
## This is DHARMa 0.4.3. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa') Note: Syntax of plotResiduals has changed in 0.3.0, see ?plotResiduals for details
library(Hmisc)
## Warning: package 'Hmisc' was built under R version 4.1.3
## Loading required package: survival
## Loading required package: Formula
## Warning: package 'Formula' was built under R version 4.1.1
## 
## Attaching package: 'Hmisc'
## The following objects are masked from 'package:dplyr':
## 
##     src, summarize
## The following objects are masked from 'package:base':
## 
##     format.pval, units
library(phia)
## Warning: package 'phia' was built under R version 4.1.1
library(lsmeans)
## Warning: package 'lsmeans' was built under R version 4.1.1
## Loading required package: emmeans
## Warning: package 'emmeans' was built under R version 4.1.1
## The 'lsmeans' package is now basically a front end for 'emmeans'.
## Users are encouraged to switch the rest of the way.
## See help('transition') for more information, including how to
## convert old 'lsmeans' objects and scripts to work with 'emmeans'.
library(emmeans)
library(multcomp)
## Warning: package 'multcomp' was built under R version 4.1.1
## Loading required package: mvtnorm
## Warning: package 'mvtnorm' was built under R version 4.1.1
## Loading required package: TH.data
## Warning: package 'TH.data' was built under R version 4.1.1
## Loading required package: MASS
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
## 
## Attaching package: 'TH.data'
## The following object is masked from 'package:MASS':
## 
##     geyser
library(plotly)
## Warning: package 'plotly' was built under R version 4.1.1
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:MASS':
## 
##     select
## The following object is masked from 'package:Hmisc':
## 
##     subplot
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
library(lmerTest)
## Warning: package 'lmerTest' was built under R version 4.1.1
## 
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
## 
##     lmer
## The following object is masked from 'package:stats':
## 
##     step
library(readxl)
library(tinytex)
## Warning: package 'tinytex' was built under R version 4.1.1
library(tidyverse)
## Warning: package 'tidyverse' was built under R version 4.1.1
## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --
## v tibble  3.1.6     v purrr   0.3.4
## v tidyr   1.2.0     v stringr 1.4.0
## v readr   2.0.1
## Warning: package 'tibble' was built under R version 4.1.2
## Warning: package 'tidyr' was built under R version 4.1.2
## Warning: package 'readr' was built under R version 4.1.1
## Warning: package 'purrr' was built under R version 4.1.1
## Warning: package 'stringr' was built under R version 4.1.1
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x tidyr::expand()    masks Matrix::expand()
## x plotly::filter()   masks dplyr::filter(), stats::filter()
## x dplyr::lag()       masks stats::lag()
## x tidyr::pack()      masks Matrix::pack()
## x dplyr::recode()    masks car::recode()
## x plotly::select()   masks MASS::select(), dplyr::select()
## x purrr::some()      masks car::some()
## x Hmisc::src()       masks dplyr::src()
## x Hmisc::summarize() masks dplyr::summarize()
## x tidyr::unpack()    masks Matrix::unpack()
library(devtools)
## Warning: package 'devtools' was built under R version 4.1.2
## Loading required package: usethis
## Warning: package 'usethis' was built under R version 4.1.2
## 
## Attaching package: 'devtools'
## The following object is masked from 'package:emmeans':
## 
##     test
library(brms)
## Warning: package 'brms' was built under R version 4.1.2
## Loading required package: Rcpp
## Warning: package 'Rcpp' was built under R version 4.1.1
## Loading 'brms' package (version 2.16.3). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').
## 
## Attaching package: 'brms'
## The following object is masked from 'package:survival':
## 
##     kidney
## The following object is masked from 'package:lme4':
## 
##     ngrps
## The following object is masked from 'package:stats':
## 
##     ar
library(ggthemes)
## Warning: package 'ggthemes' was built under R version 4.1.2

###Set working directory

#Original
#setwd("C:/Users/riant/Documents/Bachelor thesis/Data")

#Russell's PC
setwd("C:/Users/RussellChan/OneDrive - University of Twente/2020_MSCA_IF/2_Bachelor_thesis/Topics_2022/Dataframe/For Analysis/Rian Analysis")

##Import dataset

df_FAMOM <- read_xlsx("FAMOMCON_Final.xlsx")

###Filtering data
df_FAM <- filter(df_FAMOM, Group %in% c("FAM","CON"))

df_FAM_acc <- read_xlsx("sm acc_1.xlsx")

df_FAM_acc <- df_FAM_acc %>% select(-(Row_Labels:Sum_8))

df_FAM_all <- merge(df_FAM_acc, df_FAM, by = c("Subject","Session"))


df_FAM_testing <- subset(df_FAM_all, Session > 6)
df_FAM_training <- subset(df_FAM_all, Session < 7)

###############For ACC analysis######################
df_ACCfinal <- read.table("ACCoverall.csv", sep = ",", header = T, stringsAsFactors = F)

df_FAM.ACC_testing <- subset(df_ACCfinal, Session > 6)
df_FAM.ACC_training <- subset(df_ACCfinal, Session < 7)
###Creating factors for ACC
##For training
df_FAM.ACC_training$Subject <- factor(df_FAM.ACC_training$Subject)
df_FAM.ACC_training$Group <- factor(df_FAM.ACC_training$Group)
df_FAM.ACC_training$Session <- as.numeric(df_FAM.ACC_training$Session, ordered = TRUE, levels=c('1', '2', '3', '4', '5', '6'))
#df_FAM.ACC_training$Session <- factor(df_FAM.ACC_training$Session, ordered = TRUE, levels=c('1', '2', '3', '4', '5', '6'))

###Accuracy

#Training
m.df_training.ACC <- lmer(Percentage ~ Group * Session +  (1|Subject), data=df_FAM.ACC_training, REML = FALSE)
Anova(m.df_training.ACC)
summary(m.df_training.ACC, ddf="Satterthwaite")
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Percentage ~ Group * Session + (1 | Subject)
##    Data: df_FAM.ACC_training
## 
##      AIC      BIC   logLik deviance df.resid 
##    917.1    935.0   -452.6    905.1      138 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.7062 -0.4714  0.1647  0.6717  1.8075 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Subject  (Intercept)  7.262   2.695   
##  Residual             26.755   5.173   
## Number of obs: 144, groups:  Subject, 24
## 
## Fixed effects:
##                  Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)       88.9062     1.5929 107.2228  55.813   <2e-16 ***
## GroupFAM          -1.4043     2.2527 107.2228  -0.623   0.5344    
## Session            0.8342     0.3569 120.0000   2.337   0.0211 *  
## GroupFAM:Session   0.6300     0.5048 120.0000   1.248   0.2145    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) GrpFAM Sessin
## GroupFAM    -0.707              
## Session     -0.784  0.555       
## GrpFAM:Sssn  0.555 -0.784 -0.707
####estimated marginal means for significant interactions <- NOT doing it.
#emmeans(m.df_training.ACC, pairwise ~ Group | Session)

ae.m.training.ACC <- allEffects(m.df_training.ACC)
ae.m.df.training.ACC <- as.data.frame(ae.m.training.ACC[[1]])
ae.m.df.training.ACC$Session <- as.numeric(ae.m.df.training.ACC$Session, ordered = TRUE, levels=c('1', '2', '3', '4', '5', '6'))

Training.ACC <- ggplot(ae.m.df.training.ACC, aes(x=Session, y=fit, ymin=fit-se, ymax=fit+se)) + 
  geom_point(aes(color = Group)) +
  geom_path(aes(x=Session, y=fit, group=Group, colour=Group)) +
  geom_ribbon(aes(ymin=fit-se, ymax=fit+se, group=Group, fill=Group), alpha=0.2) +
  ylab("Accuracy (%)")+
  scale_x_continuous(name="Session", limits = c(1, 6), breaks=c(1,2,3,4,5,6))+
  ggtitle("Accuracy (%): Session x Group")+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))

print(Training.ACC)

###Filtering out inaccurate trials
df_FAM_acc2 <- df_FAM_all %>% filter(feedback.ACC.trial==1)


##Finding the mean and SD for each block

#1

block_1 <-filter(df_FAM_acc2, Session == 1)

mean(block_1$feedback.RT.mean) 
## [1] 620.335
sd(block_1$feedback.RT.mean)
## [1] 337.8138
#2

block_2 <-filter(df_FAM_acc2, Session == 2)

mean(block_2$feedback.RT.mean) 
## [1] 418.8285
sd(block_2$feedback.RT.mean)
## [1] 199.8392
#3

block_3 <-filter(df_FAM_acc2, Session == 3)

mean(block_3$feedback.RT.mean) 
## [1] 386.1655
sd(block_3$feedback.RT.mean)
## [1] 213.2458
#4

block_4 <-filter(df_FAM_acc, Session == 4)

mean(block_4$feedback.RT.mean) 
## Warning: Unknown or uninitialised column: `feedback.RT.mean`.
## Warning in mean.default(block_4$feedback.RT.mean): argument is not numeric or
## logical: returning NA
## [1] NA
sd(block_4$feedback.RT.mean)
## Warning: Unknown or uninitialised column: `feedback.RT.mean`.
## [1] NA
#5

block_5 <-filter(df_FAM_acc2, Session == 5)

mean(block_5$feedback.RT.mean) 
## [1] 356.5704
sd(block_5$feedback.RT.mean)
## [1] 169.9788
#6

block_6 <-filter(df_FAM_acc2, Session == 6)

mean(block_6$feedback.RT.mean) 
## [1] 348.0109
sd(block_6$feedback.RT.mean)
## [1] 165.1839
###Correct
#7

block_7 <-filter(df_FAM_acc2, Session == 7)

mean(block_7$feedback.RT.mean) 
## [1] 366.4313
sd(block_7$feedback.RT.mean)
## [1] 175.2055
#8

block_8 <-filter(df_FAM_acc2, Session == 8)

mean(block_8$feedback.RT.mean) 
## [1] 387.0527
sd(block_8$feedback.RT.mean)
## [1] 178.9957
###Filtering out the values 
data_1 <- filter(df_FAM_acc2, Session == 1 & feedback.RT.mean < 1464)
data_2 <- filter(df_FAM_acc2, Session == 2 & feedback.RT.mean < 918)
data_3 <- filter(df_FAM_acc2, Session == 3 & feedback.RT.mean < 919)
data_4 <- filter(df_FAM_acc2, Session == 4 & feedback.RT.mean < 795)
data_5 <- filter(df_FAM_acc2, Session == 5 & feedback.RT.mean < 781)
data_6 <- filter(df_FAM_acc2, Session == 6 & feedback.RT.mean < 760)

#Subset earlier and then base it on FAM/UNFAM
############
#data_7 <- filter(df_FAM_acc2, Sequence == FAM & feedback.RT.mean < 804)
#data_8 <- filter(df_FAM_acc2, Sequence == UNFAM & feedback.RT.mean < 834)


df_FAM2_training <- rbind(data_1, data_2, data_3, data_4, data_5, data_6)
#df_FAM2 <- rbind(data_1, data_2, data_3, data_4, data_5, data_6, data_7, data_8)

###Filtering training and testing
#df_FAM2_testing <- subset(df_FAM2, Session > 6)
#df_FAM2_training <- subset(df_FAM2, Session < 7)

###Creating new factors

##For all data
#df_FAM2$Subject <- factor(df_FAM2$Subject)
#df_FAM2$Group <- factor(df_FAM2$Group)
#df_FAM2$FAM_UNFAM <- factor(df_FAM2$FAM_UNFAM)

##For training
df_FAM2_training$Subject <- factor(df_FAM2_training$Subject)
df_FAM2_training$Group <- factor(df_FAM2_training$Group)
df_FAM2_training$Session <- as.numeric(df_FAM2_training$Session, ordered = TRUE, levels=c('1', '2', '3', '4', '5', '6'))
#df_FAM2_training$Session <- factor(df_FAM2_training$Session, ordered = TRUE, levels=c('1', '2', '3', '4', '5', '6'))

###Mean RT

##Training
m.df_training.meanRT <- lmer(feedback.RT.mean ~ Group * Session +  (1|Subject), data=df_FAM2_training, REML = FALSE)
Anova(m.df_training.meanRT)
summary(m.df_training.meanRT, ddf="Satterthwaite")
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: feedback.RT.mean ~ Group * Session + (1 | Subject)
##    Data: df_FAM2_training
## 
##      AIC      BIC   logLik deviance df.resid 
##  72642.3  72682.2 -36315.1  72630.3     5775 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.7050 -0.5648 -0.1305  0.3573  7.3528 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Subject  (Intercept) 20743    144.0   
##  Residual             16350    127.9   
## Number of obs: 5781, groups:  Subject, 24
## 
## Fixed effects:
##                  Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)       517.026     41.954   24.708  12.324 4.78e-12 ***
## GroupFAM           37.534     59.339   24.719   0.633    0.533    
## Session           -33.491      1.430 5757.105 -23.427  < 2e-16 ***
## GroupFAM:Session   -8.955      2.023 5757.128  -4.428 9.71e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) GrpFAM Sessin
## GroupFAM    -0.707              
## Session     -0.121  0.086       
## GrpFAM:Sssn  0.086 -0.122 -0.707
####estimated marginal means for significant interactions
#emmeans(m.df_training.meanRT, pairwise ~ Group | Session)

ae.m.training.meanRT <- allEffects(m.df_training.meanRT)
ae.m.df.training.meanRT <- as.data.frame(ae.m.training.meanRT[[1]])
#ae.m.df.training.meanRT$Session <- factor(ae.m.df.training.meanRT$Session, ordered = TRUE, levels=c('1', '2', '3', '4', '5', '6'))

#OLD plot.  Still works
Training.meanRT <- ggplot(ae.m.df.training.meanRT, aes(x=Session, y=fit, color=Group)) + 
  geom_errorbar(aes(ymin=lower, ymax=upper), width=.1) +
  geom_line() +
  geom_point() +
  ylab("Mean reaction time (ms)") +
  xlab("Session") +
  ggtitle("Mean reaction time (ms) by Group") +
  theme_classic()

print(Training.meanRT)

#Here we used standard error of the mean for the plot.  It looks nicer.
ae <- ggplot(ae.m.df.training.meanRT, aes(x=Session, y=fit, ymin=lower, ymax=upper))+
  geom_point(aes(color = Group)) +
  geom_path(aes(x=Session, y=fit, group=Group, colour=Group)) +
  geom_ribbon(aes(ymin=fit-se, ymax=fit+se, group=Group, fill=Group), alpha=0.2) +
  ylab("Mean reaction time (ms)")+
  scale_x_continuous(name="Session", limits = c(1, 6), breaks=c(1,2,3,4,5,6))+
  ggtitle("Mean reaction time (ms): Session x Group")+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))
plot(ae)

#Explanation
#A linear mixed effects model was run on RT (ms), using a fixed effect for the numeric predictor variable of training sessions and random effects for the subject variable.  The results reveled no main effect of Group (p = .93; FAM or Control), a significant main linear effect of Session (Block 1 to 6) on RT, χ²(1, N = 24) = 1409.4, p<.001, reflecting shorter RT across the training blocks for both groups. In addition, a significant Group x Session interaction, χ²(1, N = 24) = 19.6, p<.001.  This model show than the linear slopes between the FAM and Control Group were significantly different and that the FAM group were predicted to have shorter RTs with training compared to the Control group.

#For 1st part of the discussion:
#In this case, the null hypothesis of FAM and Control group progressing similarly in learning phase of sequence learning can be rejected and conclude that FAM