##Set working directory

setwd("C:/Users/ChanRW/OneDrive - Universiteit Twente/2020_Alpha_SeqLearn/Current Analysis")

##Packages

library(lme4)
library(effects)
library(lattice)
library(car)
library(ggplot2) 
library(knitr)
library(reshape2)
library(dplyr)
library(forcats)
library(DHARMa)
library(Hmisc)
#library(qpcR)
library(phia)
library(lsmeans)
library(emmeans)
library(multcomp)
library(plotly)
library(lmerTest)

##Import dataset

d.MU<-read.table("Pow_Data_DEC18_V3.4.CSV", sep = ",", header = T, stringsAsFactors = F)

#Factors

#Create Factors
d.MU$ID <- factor(d.MU$ID)
d.MU$Group <- factor(d.MU$Group)
d.MU$RoI <- factor(d.MU$RoI, levels=c('F', 'C', 'P'))
d.MU$Session <- factor(d.MU$Session, levels=c('1', '2', '3'))
d.MU$Timepoint <-as.numeric(d.MU$Timepoint, levels=c('1', '2', '3', '4', '5', '6', '7', '8', '9', '10'))
#d.MU$Timepoint <- factor(d.MU$Timepoint)
d.MU$muTheta <- as.numeric(d.MU$muTheta)
d.MU$muAlpha <- as.numeric(d.MU$muAlpha)

#Model1: Theta Random Intercepts: ID and session(timepoint). More conservative, but this might be the most accurate for our design since we are interested in Sessions and Timepoints muTheta changes as an effect of the Group. Timepoints has to be a numeric for this to work.

m.T1 <- lmer(muTheta ~ Group * Session * RoI  + (1 | ID) + (1 | Session:Timepoint), data=d.MU, REML = FALSE)
Anova(m.T1)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: muTheta
##                       Chisq Df Pr(>Chisq)    
## Group               12.8539  2   0.001617 ** 
## Session              7.4757  2   0.023805 *  
## RoI               7326.0047  2  < 2.2e-16 ***
## Group:Session      534.0749  4  < 2.2e-16 ***
## Group:RoI          162.9969  4  < 2.2e-16 ***
## Session:RoI         37.8648  4  1.195e-07 ***
## Group:Session:RoI  102.2639  8  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m.T1test <- as(m.T1,"lmerModLmerTest")
summary(m.T1test, ddf="Satterthwaite")
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: muTheta ~ Group * Session * RoI + (1 | ID) + (1 | Session:Timepoint)
##    Data: d.MU
## 
##       AIC       BIC    logLik  deviance  df.resid 
##  213316.0  213580.4 -106628.0  213256.0     49720 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.0963 -0.6237 -0.1151  0.4834 10.8887 
## 
## Random effects:
##  Groups            Name        Variance Std.Dev.
##  Session:Timepoint (Intercept) 0.06082  0.2466  
##  ID                (Intercept) 0.77156  0.8784  
##  Residual                      4.23529  2.0580  
## Number of obs: 49750, groups:  Session:Timepoint, 30; ID, 29
## 
## Fixed effects:
##                            Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)               5.181e+00  3.239e-01  3.395e+01  15.995  < 2e-16 ***
## GroupMED1                 7.806e-01  4.321e-01  3.027e+01   1.807 0.080789 .  
## GroupMED21                8.233e-01  4.059e-01  3.027e+01   2.028 0.051401 .  
## Session2                  1.653e-01  1.303e-01  5.230e+01   1.268 0.210308    
## Session3                  6.333e-02  1.303e-01  5.230e+01   0.486 0.628994    
## RoIC                     -1.177e+00  7.327e-02  4.969e+04 -16.069  < 2e-16 ***
## RoIP                     -2.343e+00  7.356e-02  4.969e+04 -31.846  < 2e-16 ***
## GroupMED1:Session2        1.793e-01  9.535e-02  4.969e+04   1.881 0.060000 .  
## GroupMED21:Session2       3.300e-01  9.056e-02  4.969e+04   3.644 0.000269 ***
## GroupMED1:Session3       -7.338e-02  9.535e-02  4.969e+04  -0.770 0.441532    
## GroupMED21:Session3       6.227e-01  8.963e-02  4.969e+04   6.947 3.78e-12 ***
## GroupMED1:RoIC            6.358e-01  1.006e-01  4.969e+04   6.320 2.64e-10 ***
## GroupMED21:RoIC           2.755e-01  9.455e-02  4.969e+04   2.914 0.003570 ** 
## GroupMED1:RoIP            8.744e-01  1.010e-01  4.969e+04   8.659  < 2e-16 ***
## GroupMED21:RoIP           3.595e-01  9.482e-02  4.969e+04   3.791 0.000150 ***
## Session2:RoIC            -2.864e-01  1.036e-01  4.969e+04  -2.764 0.005704 ** 
## Session3:RoIC            -1.338e-01  1.035e-01  4.969e+04  -1.292 0.196196    
## Session2:RoIP            -6.035e-02  1.038e-01  4.969e+04  -0.581 0.561058    
## Session3:RoIP             1.789e-01  1.037e-01  4.969e+04   1.725 0.084554 .  
## GroupMED1:Session2:RoIC   8.258e-03  1.423e-01  4.969e+04   0.058 0.953737    
## GroupMED21:Session2:RoIC  9.178e-02  1.349e-01  4.969e+04   0.680 0.496289    
## GroupMED1:Session3:RoIC  -2.925e-01  1.423e-01  4.969e+04  -2.056 0.039798 *  
## GroupMED21:Session3:RoIC -1.800e-02  1.336e-01  4.969e+04  -0.135 0.892852    
## GroupMED1:Session2:RoIP   1.504e-01  1.425e-01  4.969e+04   1.055 0.291299    
## GroupMED21:Session2:RoIP  1.317e-02  1.351e-01  4.969e+04   0.097 0.922350    
## GroupMED1:Session3:RoIP  -7.264e-01  1.425e-01  4.969e+04  -5.098 3.44e-07 ***
## GroupMED21:Session3:RoIP  2.641e-01  1.338e-01  4.969e+04   1.973 0.048451 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 27 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it

#Model2: Alpha Random Intercepts: ID and session(timepoint). More conservative, but this might be the most accurate for our design since we are interested in Sessions and Timepoints muAlpha changes as an effect of the Group. Timepoints has to be a numeric for this to work.

m.A1 <- lmer(muAlpha ~ Group * Session * RoI  + (1 | ID) + (1 | Session:Timepoint), data=d.MU, REML = FALSE)
Anova(m.A1)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: muAlpha
##                        Chisq Df Pr(>Chisq)    
## Group                 6.6552  2    0.03588 *  
## Session              18.4274  2  9.967e-05 ***
## RoI               45626.2874  2  < 2.2e-16 ***
## Group:Session       159.8570  4  < 2.2e-16 ***
## Group:RoI          2397.3415  4  < 2.2e-16 ***
## Session:RoI           2.6630  4    0.61571    
## Group:Session:RoI   174.9800  8  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m.A1test <- as(m.A1,"lmerModLmerTest")
summary(m.A1test, ddf="Satterthwaite")
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: muAlpha ~ Group * Session * RoI + (1 | ID) + (1 | Session:Timepoint)
##    Data: d.MU
## 
##       AIC       BIC    logLik  deviance  df.resid 
##  247728.2  247992.6 -123834.1  247668.2     49720 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.2652 -0.6447 -0.0277  0.6345  5.1216 
## 
## Random effects:
##  Groups            Name        Variance Std.Dev.
##  Session:Timepoint (Intercept) 0.07988  0.2826  
##  ID                (Intercept) 4.81715  2.1948  
##  Residual                      8.45472  2.9077  
## Number of obs: 49750, groups:  Session:Timepoint, 30; ID, 29
## 
## Fixed effects:
##                            Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)                  5.8555     0.7842    30.1759   7.467 2.44e-08 ***
## GroupMED1                    0.1441     1.0707    29.4034   0.135  0.89389    
## GroupMED21                  -0.5418     1.0058    29.4041  -0.539  0.59413    
## Session2                     0.8586     0.1599    65.8039   5.368 1.11e-06 ***
## Session3                     0.8116     0.1599    65.8039   5.074 3.41e-06 ***
## RoIC                         5.3428     0.1035 49692.0126  51.611  < 2e-16 ***
## RoIP                         9.1458     0.1039 49692.0291  87.997  < 2e-16 ***
## GroupMED1:Session2          -0.9832     0.1347 49692.0105  -7.299 2.95e-13 ***
## GroupMED21:Session2         -0.2480     0.1280 49692.5190  -1.938  0.05262 .  
## GroupMED1:Session3          -1.1024     0.1347 49692.0106  -8.183 2.82e-16 ***
## GroupMED21:Session3          0.1898     0.1266 49692.0106   1.499  0.13396    
## GroupMED1:RoIC              -3.2890     0.1421 49692.0116 -23.139  < 2e-16 ***
## GroupMED21:RoIC             -2.8276     0.1336 49692.0121 -21.167  < 2e-16 ***
## GroupMED1:RoIP              -3.3668     0.1427 49692.0218 -23.597  < 2e-16 ***
## GroupMED21:RoIP             -3.1709     0.1340 49692.0217 -23.668  < 2e-16 ***
## Session2:RoIC               -0.3975     0.1464 49692.0128  -2.715  0.00662 ** 
## Session3:RoIC                0.1259     0.1463 49692.0116   0.861  0.38924    
## Session2:RoIP               -0.6935     0.1467 49692.0217  -4.727 2.28e-06 ***
## Session3:RoIP               -0.2462     0.1466 49692.0199  -1.680  0.09303 .  
## GroupMED1:Session2:RoIC      0.4887     0.2011 49692.0121   2.430  0.01510 *  
## GroupMED21:Session2:RoIC     0.5944     0.1906 49692.0120   3.118  0.00182 ** 
## GroupMED1:Session3:RoIC      0.1618     0.2010 49692.0115   0.805  0.42081    
## GroupMED21:Session3:RoIC    -0.3757     0.1888 49692.0111  -1.990  0.04663 *  
## GroupMED1:Session2:RoIP      0.6373     0.2014 49692.0172   3.164  0.00156 ** 
## GroupMED21:Session2:RoIP     0.9735     0.1909 49692.0171   5.101 3.40e-07 ***
## GroupMED1:Session3:RoIP      1.2119     0.2013 49692.0162   6.020 1.75e-09 ***
## GroupMED21:Session3:RoIP    -0.4795     0.1891 49692.0164  -2.536  0.01121 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 27 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it

#Modelfit checking.

AIC(m.T1, m.A1)
##      df      AIC
## m.T1 30 213316.0
## m.A1 30 247728.2

#All effects Model Alpha Theta.

#Model 1 Effects
ae.m.T1 <- allEffects(m.T1)
ae.m.df.T1 <- as.data.frame(ae.m.T1[[1]])
#Ordering RoI
ae.m.df.T1$RoI <- factor(ae.m.df.T1$RoI, levels=c('F', 'C', 'P'))
#Ordering Timepoint
ae.m.df.T1$Session <- as.character(ae.m.df.T1$Session)
ae.m.df.T1$Session <- as.numeric(ae.m.df.T1$Session)

#Model 1 Effects
ae.m.A1 <- allEffects(m.A1)
ae.m.df.A1 <- as.data.frame(ae.m.A1[[1]])
#Ordering RoI
ae.m.df.A1$RoI <- factor(ae.m.df.A1$RoI, levels=c('F', 'C', 'P'))
#Ordering Timepoint
ae.m.df.A1$Session <- as.character(ae.m.df.A1$Session)
ae.m.df.A1$Session <- as.numeric(ae.m.df.A1$Session)

#Plotting models 0 to 4. They all look the same….

#Model T1 plot: with 95% CIs
#ae.T1 <- ggplot(ae.m.df.T1, aes(x=Session,y=fit, color=Group))+
  #geom_line() +
  #geom_path(aes(x=Session, y=fit, color=Group)) +
  #geom_errorbar(aes(ymin=lower, ymax=upper), width=.1) +
  #geom_point(aes(color = Group))+
  #ylab("Mean Theta Power (n.u.)")+
  #scale_x_continuous(name="Session", breaks=c(1,2,3) , limits=c(0.75,3.25))+
  #ggtitle("muTheta Group x Session x RoI")+
  #theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = #element_line(colour = "black"))+
  #facet_grid(RoI~.)
#plot(ae.T1)

#This with Ribbon
#Model T1 plot: with 95% CIs
ae.T1 <- ggplot(ae.m.df.T1, aes(x=Session,y=fit, group=Group))+
  geom_ribbon(aes(ymin=lower, ymax=upper, fill=Group), alpha=0.2) +
  geom_line(aes(color=Group),size=1) +
  geom_point(aes(color = Group, shape = Group), size=2.2)+
  ylab("Mean Theta Power (n.u.)")+
  scale_x_continuous(name="Session", breaks=c(1,2,3) , limits=c(0.75,3.25))+
  ggtitle("muTheta Group x Session x RoI")+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"), axis.text = element_text(size = 25))+
  facet_grid(RoI~.)
plot(ae.T1)

#Copied from another plot
#ae.position<-ggplot(ae.m.df.1, aes(x=Time,y=fit, group=Block))+
  #geom_ribbon(aes(ymin=lower, ymax=upper, fill=Block), alpha=0.2) +
  #geom_line(aes(size=0.5, color=Block)) +
  #geom_point(aes(color=Block, size=2))+
  #ylab("Relative β ERDS Power (%)")+
  #xlab("Timepoints")+
  #theme_classic()

#Printing Session effects facet
#print(ae.position)


#Model T1.1 plot: With fit+se
ae.T1.1 <- ggplot(ae.m.df.T1, aes(x=Session,y=fit, color=Group))+
  geom_line() +
  geom_path(aes(x=Session, y=fit, color=Group)) +
  geom_errorbar(aes(ymin=fit-se, ymax=fit+se), width=.1) +
  geom_point(aes(color = Group))+
  ylab("Mean Theta Power (n.u.)")+
  scale_x_continuous(name="Session", breaks=c(1,2,3) , limits=c(0.75,3.25))+
  ggtitle("Theta Group x Session x RoI")+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"), axis.text = element_text(size = 20))+
  facet_grid(RoI~.)
plot(ae.T1.1)

#Model A1 plot: with 95% CIs
#ae.A1 <- ggplot(ae.m.df.A1, aes(x=Session,y=fit, color=Group))+
  #geom_line() +
  #geom_path(aes(x=Session, y=fit, color=Group)) +
  #geom_errorbar(aes(ymin=lower, ymax=upper), width=.1) +
  #geom_point(aes(color = Group))+
  #ylab("Mean Theta Power (n.u.)")+
  #scale_x_continuous(name="Session", breaks=c(1,2,3) , limits=c(0.75,3.25))+
  #ggtitle("muTheta Group x Session x RoI")+
  #theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))+
  #facet_grid(RoI~.)
#plot(ae.A1)

#This with Ribbon
#Model T1 plot: with 95% CIs
ae.A1 <- ggplot(ae.m.df.A1, aes(x=Session,y=fit, group=Group))+
  geom_ribbon(aes(ymin=lower, ymax=upper, fill=Group), alpha=0.2) +
  geom_line(aes(color=Group),size=1) +
  geom_point(aes(color = Group, shape = Group), size=2.2)+
  ylab("Mean Alpha Power (n.u.)")+
  scale_x_continuous(name="Session", breaks=c(1,2,3) , limits=c(0.75,3.25))+
  ggtitle("Alpha Group x Session x RoI")+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"), axis.text = element_text(size = 25))+
  facet_grid(RoI~.)
plot(ae.A1)

#Model 1.1 plot: With fit+se
ae.A1.1 <- ggplot(ae.m.df.A1, aes(x=Session,y=fit, color=Group))+
  geom_line() +
  geom_path(aes(x=Session, y=fit, color=Group)) +
  geom_errorbar(aes(ymin=fit-se, ymax=fit+se), width=.1) +
  geom_point(aes(color = Group))+
  ylab("Mean Theta Power (n.u.)")+
  scale_x_continuous(name="Session", breaks=c(1,2,3) , limits=c(0.75,3.25))+
  ggtitle("muTheta Group x Session x RoI")+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))+
  facet_grid(RoI~.)
plot(ae.A1.1)

print(ae.T1) print(ae.T1.1) print(ae.A1) print(ae.A1.1)

####Posthocs? I don’t think we need it as we should rely on plots.###

#Method 1: To get summary of posthocs, simple but not the most correct since they are only conparing against Control. #summary(glht(m.ThetaRanSS1, linfct = mcp(Group = “Tukey”)), test = adjusted(“holm”)) #lsmeans(m.ThetaRanSS1, pairwise~Group, adjust=“tukey”) #emmeans(m.ThetaRanSS1, pairwise~Group)

#summary(glht(m.ThetaRanSS2, linfct = mcp(Group = “Tukey”)), test = adjusted(“holm”)) #lsmeans(m.ThetaRanSS2, pairwise~GroupSession, adjust=“tukey”) #emmeans(m.ThetaRanSS2, pairwise~GroupSession)

#lsmeans(m.ThetaRanSS3, pairwise~GroupSessionTimepoint, adjust=“tukey”) #emmeans(m.ThetaRanSS3, pairwise~GroupSessionTimepoint)

#lsmeans(m.ThetaRanSS4, pairwise~GroupSessionTimepointRoI, adjust=“tukey”) #emmeans(m.ThetaRanSS3, pairwise~GroupSession*Timepoint)

#Method 2: To get summary of posthocs - This is the more correct way. #emmeans(MODEL_NAME, pairwise ~ size) #Overall emmeans -> Hasty and not correct #Create new interaction matrices #EG. emm_NAME <- emmeans(MODEL_NAME, pairwise ~ size | type) #The last term is the conditioning on type.