Broncho_Behav V1.0
Russell Chan
21 Dec 2021
##Set working directory
#Home
setwd("C:/Users/RussellChan/OneDrive - University of Twente/2020_MSCA_IF/2_Bachelor_thesis/2022_Joel_Broncho_X/Data")
#Work
#setwd("C:/Users/ChanRW/OneDrive - Universiteit Twente/2021_Oscillatory EXP 3 Paper 4 Analysis/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(phia)
library(lsmeans)
library(emmeans)
library(multcomp)
library(plotly)
library(lmerTest)
library(readxl)
library(tinytex)##Import dataset
df1<-read.table("Task1Behav.csv", sep = ",", header = T, stringsAsFactors = F)
df3<-read.table("Task3Behav.csv", sep = ",", header = T, stringsAsFactors = F)#Factors
#Create Factors
df1$Participant <- factor(df1$Participant)
df1$Trial <- factor(df1$Trial, ordered = TRUE, levels=c('1', '2', '3', '4', '5'))
df3$Participant <- factor(df3$Participant)
df3$Trial <- factor(df3$Trial, ordered = TRUE, levels=c('1', '2', '3'))#Task 1 Model
#Task 1 Model
m.df1 <- lmer(Movement_Time ~ Trial + (1|Participant), data=df1, REML = FALSE)
Anova(m.df1)summary(m.df1, ddf="Satterthwaite")## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: Movement_Time ~ Trial + (1 | Participant)
## Data: df1
##
## AIC BIC logLik deviance df.resid
## 532.4 545.8 -259.2 518.4 43
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.7814 -0.5893 -0.1268 0.3660 2.3728
##
## Random effects:
## Groups Name Variance Std.Dev.
## Participant (Intercept) 747.2 27.33
## Residual 1442.0 37.97
## Number of obs: 50, groups: Participant, 10
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 102.034 10.176 10.000 10.026 1.55e-06 ***
## Trial.L -48.326 12.009 40.000 -4.024 0.000247 ***
## Trial.Q 24.501 12.009 40.000 2.040 0.047964 *
## Trial.C -34.722 12.009 40.000 -2.891 0.006172 **
## Trial^4 7.499 12.009 40.000 0.624 0.535883
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Tril.L Tril.Q Tril.C
## Trial.L 0.000
## Trial.Q 0.000 0.000
## Trial.C 0.000 0.000 0.000
## Trial^4 0.000 0.000 0.000 0.000######################
#Task 3 Model
m.df3 <- lmer(Movement_Time ~ Trial + (1|Participant), data=df3, REML = FALSE)
Anova(m.df3)summary(m.df3, ddf="Satterthwaite")## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: Movement_Time ~ Trial + (1 | Participant)
## Data: df3
##
## AIC BIC logLik deviance df.resid
## 364.5 371.5 -177.3 354.5 25
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.58011 -0.63303 -0.01445 0.45310 2.34870
##
## Random effects:
## Groups Name Variance Std.Dev.
## Participant (Intercept) 12256 110.70
## Residual 3522 59.35
## Number of obs: 30, groups: Participant, 10
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 533.52 36.65 10.00 14.559 4.66e-08 ***
## Trial.L -116.71 18.77 20.00 -6.219 4.49e-06 ***
## Trial.Q 24.39 18.77 20.00 1.299 0.209
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Tril.L
## Trial.L 0.000
## Trial.Q 0.000 0.000plot(m.df1)
#Alpha Posthocs to determine exact locus
#Some options
#Holm Contrasts
#https://stats.stackexchange.com/questions/237512/how-to-perform-post-hoc-test-on-lmer-model
#Interpretation from Summary of T values
#https://stats.stackexchange.com/questions/87412/how-to-interpret-2-way-and-3-way-interaction-in-lmer/87415
#summary(glht(m.AF_WO, linfct = mcp(Group = "Tukey")), test = adjusted("holm"))
#plot(emmeans(m.AF_WO2, ~ Group | Session))
#################################
#Task 1 model post hoc
summary(glht(m.df1)) ##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: lmer(formula = Movement_Time ~ Trial + (1 | Participant), data = df1,
## REML = FALSE)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) == 0 102.034 10.176 10.026 < 2e-16 ***
## Trial.L == 0 -48.326 12.009 -4.024 0.000286 ***
## Trial.Q == 0 24.501 12.009 2.040 0.190232
## Trial.C == 0 -34.722 12.009 -2.891 0.019025 *
## Trial^4 == 0 7.499 12.009 0.624 0.977630
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)emmeans(m.df1, pairwise ~ Trial)## $emmeans
## Trial emmean SE df lower.CL upper.CL
## 1 157.6 15.6 37.9 126.0 189
## 2 85.2 15.6 37.9 53.6 117
## 3 94.3 15.6 37.9 62.7 126
## 4 98.6 15.6 37.9 67.0 130
## 5 74.5 15.6 37.9 42.9 106
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## 1 - 2 72.35 17.9 44.4 4.042 0.0019
## 1 - 3 63.26 17.9 44.4 3.534 0.0082
## 1 - 4 58.99 17.9 44.4 3.295 0.0158
## 1 - 5 83.09 17.9 44.4 4.641 0.0003
## 2 - 3 -9.09 17.9 44.4 -0.508 0.9862
## 2 - 4 -13.36 17.9 44.4 -0.746 0.9443
## 2 - 5 10.74 17.9 44.4 0.600 0.9744
## 3 - 4 -4.26 17.9 44.4 -0.238 0.9993
## 3 - 5 19.83 17.9 44.4 1.108 0.8015
## 4 - 5 24.10 17.9 44.4 1.346 0.6644
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 5 estimates#Task 3 model post hoc
summary(glht(m.df3)) ##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: lmer(formula = Movement_Time ~ Trial + (1 | Participant), data = df3,
## REML = FALSE)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) == 0 533.53 36.65 14.559 < 2e-16 ***
## Trial.L == 0 -116.71 18.77 -6.219 1.5e-09 ***
## Trial.Q == 0 24.39 18.77 1.299 0.476
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)emmeans(m.df3, pairwise ~ Trial)## $emmeans
## Trial emmean SE df lower.CL upper.CL
## 1 626 41.9 15.1 537 715
## 2 514 41.9 15.1 424 603
## 3 461 41.9 15.1 372 550
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## 1 - 2 112.4 28 22.2 4.017 0.0016
## 1 - 3 165.1 28 22.2 5.900 <.0001
## 2 - 3 52.7 28 22.2 1.882 0.1672
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimates#Alpha models effects plot
#Options
#Ordering timepoint if needed
#ae.m.df.AF2$Time <- factor(ae.m.df.AF2$Time, levels=c('1', '2', '3', '4', '5', '6', '7', '8', '9', '10'))
#Task 1 Model effects
ae.m.df1 <- allEffects(m.df1)
ae.m.df.df1 <- as.data.frame(ae.m.df1[[1]])
##########################
#Task 3 Model effects
ae.m.df3 <- allEffects(m.df3)
ae.m.df.df3 <- as.data.frame(ae.m.df3[[1]])#Alpha models effects plot
#Options: Some ggplot2 customization to implement if needed.
#limits=c(0.75,3.25)
#breaks=c(1,2,3,4,5,6,7,8,9,10)
#+facet_grid(Session~.)
#geom_line(aes(color=Group),size=1)
#geom_ribbon(aes(ymin=fit-se, ymax=fit+se, fill=Group), alpha=0.2)
#geom_smooth(aes(color=Group),size=1)
#geom_line(aes(color=Group),size=1)
#Task 1 plot
ae.Task1<-ggplot(ae.m.df.df1, aes(x=Trial,y=fit,color=Trial))+
geom_errorbar(aes(ymin=lower, ymax=upper), width=.1) +
geom_line() +
geom_point()+
ylab("Movement Time (Seconds)")+
xlab("Trial")+
ggtitle("Movement Time ~ Trial")+
theme_light()
plot(ae.Task1)## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
#Task 3 plot
ae.Task3<-ggplot(ae.m.df.df3, aes(x=Trial,y=fit,color=Trial))+
geom_errorbar(aes(ymin=lower, ymax=upper), width=.1) +
geom_line() +
geom_point()+
ylab("Movement Time (Seconds)")+
xlab("Trial")+
ggtitle("Movement Time ~ Trial")+
theme_light()
plot(ae.Task3)## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?