library(readxl)
library(lme4)
## Loading required package: Matrix
library(effects)
## Loading required package: carData
## lattice theme set by effectsTheme()
## See ?effectsTheme for details.
library(car)
library(emmeans)
library(knitr)
library(reshape2)
library(dplyr)
##
## 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)
library(DHARMa)
## This is DHARMa 0.4.6. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa')
library(Hmisc)
##
## 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)
library(lsmeans)
## 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(multcomp)
## Loading required package: mvtnorm
## Loading required package: survival
## Loading required package: TH.data
## 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(nlme)
##
## Attaching package: 'nlme'
## The following object is masked from 'package:dplyr':
##
## collapse
## The following object is masked from 'package:lme4':
##
## lmList
library(lattice)
library(ggplot2)
library(lmerTest)
##
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
##
## lmer
## The following object is masked from 'package:stats':
##
## step
library(stats)
library(fractalRegression)
library(MASS)
library(car)
setwd("C:/Users/Marcel/Desktop/Bachelor Thesis/Data Analysis/Relevant R Scripts/RT Data")
rt <-read_excel("rt_data_cleaned_merged.xlsx")
rt$Block <- factor(rt$Block, ordered = TRUE, levels = c('1','2','3','4','5','6','7','8','9','10'))
rt$Subject <- as.factor(rt$Subject)
rt$Trial <- as.numeric(rt$Trial)
str(rt)
## tibble [66,240 × 8] (S3: tbl_df/tbl/data.frame)
## $ Subject : Factor w/ 23 levels "6","11","12",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Block : Ord.factor w/ 10 levels "1"<"2"<"3"<"4"<..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Step : num [1:66240] 1 2 3 4 5 6 1 2 3 4 ...
## $ Step.ACC : num [1:66240] 1 1 1 1 1 1 1 1 1 1 ...
## $ RT : num [1:66240] 646 1252 554 2384 641 ...
## $ Trial.ACC: num [1:66240] 1 1 1 1 1 1 1 1 1 1 ...
## $ Sequence : num [1:66240] 2 2 2 2 2 2 2 2 2 2 ...
## $ Trial : num [1:66240] 1 1 1 1 1 1 2 2 2 2 ...
###### mean outlier removal
# Calculate block-level mean RT and SD for each subject, based on individual steps
block_stats_mean <- rt %>%
group_by(Subject, Block) %>%
summarise(
BlockMeanRT = mean(RT, na.rm = TRUE),
BlockSDRT = sd(RT, na.rm = TRUE),
.groups = 'drop'
)
# Join block-level stats back to the original dataset
rt_with_block_stats_mean <- rt %>%
left_join(block_stats_mean, by = c("Subject", "Block"))
# Flag steps as outliers based on the mean and 2.5*SD criterion.
rt_with_block_stats_mean <- rt_with_block_stats_mean %>%
mutate(Outlier = abs(RT - BlockMeanRT) > 2.5 * BlockSDRT)
# Create a dataset with only the individual outlier steps removed
rt_cleaned_mean_steps <- rt_with_block_stats_mean %>%
filter(!Outlier) %>%
dplyr::select(-c(BlockMeanRT, BlockSDRT, Outlier))
# Splitting the dataset into the learning blocks and test blocks
rt_learningblocks <- rt_cleaned_mean_steps %>%
filter(Block %in% c("1", "2", "3", "4", "5", "6", "7", "8"))
rt_testblocks <- rt_cleaned_mean_steps %>%
filter(Block %in% c("9", "10"))
# Calculate the mean step accuracy and mean for Response times for each block for each participant
blockaverage <- rt_learningblocks %>%
group_by(Subject, Block) %>%
summarise(MeanStepACC = mean(Step.ACC, na.rm = TRUE),
MeanRT = mean(RT, na.rm = TRUE)) %>%
ungroup()
## `summarise()` has grouped output by 'Subject'. You can override using the
## `.groups` argument.
# Plot for mean response times
rt_plot <- ggplot(blockaverage, aes(x = Block, y = MeanRT, group = Subject)) +
geom_line(aes(color = as.factor(Subject))) +
facet_wrap(~ Subject, scales = "free_x") + # Allow free scales on the x-axis only
labs(title = "A) Participants' Mean Response Time Over Blocks",
x = "Block",
y = "Mean Response Time (ms)") +
theme_minimal() +
theme(legend.position = "none") +
ylim(50, 1400) # Set fixed y-axis limits
print(rt_plot)
# Plot for mean step accuracy
accuracy_plot <- ggplot(blockaverage, aes(x = Block, y = MeanStepACC, group = Subject)) +
geom_line(aes(color = as.factor(Subject))) +
facet_wrap(~ Subject, scales = "free_x") + # Allow free scales on the x-axis only
labs(title = "B) Participants' Mean Step Accuracy Over Blocks",
x = "Block",
y = "Predicted Mean Step Accuracy") +
theme_minimal() +
theme(legend.position = "none") +
ylim(0.7, 1) # Set fixed y-axis limits
print(accuracy_plot)
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).
# Fit the linear mixed-effects model for raw response times
model_rt_lmm <- lmer(RT ~ Block + (1 | Subject), data = rt, REML = FALSE)
summary(model_rt_lmm)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: RT ~ Block + (1 | Subject)
## Data: rt
##
## AIC BIC logLik deviance df.resid
## 1052083.5 1052192.7 -526029.8 1052059.5 66228
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.493 -0.301 -0.144 0.102 108.083
##
## Random effects:
## Groups Name Variance Std.Dev.
## Subject (Intercept) 25231 158.8
## Residual 461808 679.6
## Number of obs: 66240, groups: Subject, 23
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 472.666 33.226 23.000 14.226 6.91e-13 ***
## Block.L -251.906 8.350 66217.000 -30.169 < 2e-16 ***
## Block.Q 220.624 8.350 66217.000 26.423 < 2e-16 ***
## Block.C -109.154 8.350 66217.000 -13.073 < 2e-16 ***
## Block^4 138.100 8.350 66217.000 16.540 < 2e-16 ***
## Block^5 -36.720 8.350 66217.000 -4.398 1.10e-05 ***
## Block^6 37.519 8.350 66217.000 4.494 7.02e-06 ***
## Block^7 -17.512 8.350 66217.000 -2.097 0.036 *
## Block^8 11.007 8.350 66217.000 1.318 0.187
## Block^9 2.675 8.350 66217.000 0.320 0.749
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Blck.L Blck.Q Blck.C Blck^4 Blck^5 Blck^6 Blck^7 Blck^8
## Block.L 0.000
## Block.Q 0.000 0.000
## Block.C 0.000 0.000 0.000
## Block^4 0.000 0.000 0.000 0.000
## Block^5 0.000 0.000 0.000 0.000 0.000
## Block^6 0.000 0.000 0.000 0.000 0.000 0.000
## Block^7 0.000 0.000 0.000 0.000 0.000 0.000 0.000
## Block^8 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
## Block^9 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
# Obtain effects of the RT model
effects_rt <- allEffects(model_rt_lmm)
df_rt <- as.data.frame(effects_rt[[1]])
df_rt$l83 <- df_rt$fit - 1.3722 * df_rt$se
df_rt$u83 <- df_rt$fit + 1.3722 * df_rt$se
base_size <- 14
ggplot(df_rt, aes(x = Block, y = fit)) +
geom_errorbar(aes(ymin = l83, ymax = u83), width = 0.2, size = 1, color = "darkorchid4") +
geom_point(color = "darkorchid2", size = 4) +
geom_vline(xintercept = 8.5, color = "red", size = 1) +
xlab("Block") + ylab("RT(ms)") +
ggtitle("A) RT ~ Block") +
theme_minimal() +
theme(text = element_text(size = base_size),
axis.title = element_text(size = base_size + 2),
plot.title = element_text(size = base_size + 4, face = "bold"),
axis.text = element_text(size = base_size - 2),
axis.text.y = element_text(face = "bold"),
axis.text.x = element_text(face = "bold"))
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
# Fit a logistic generalized mixed-effects model for binary step accuracy
model_acc <- glmer(Step.ACC ~ Block + (1 | Subject), data = rt, family = binomial(link = "logit"))
summary(model_acc)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Step.ACC ~ Block + (1 | Subject)
## Data: rt
##
## AIC BIC logLik deviance df.resid
## 27399.7 27499.8 -13688.8 27377.7 66229
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.7604 0.1750 0.2118 0.2560 0.7081
##
## Random effects:
## Groups Name Variance Std.Dev.
## Subject (Intercept) 0.2215 0.4706
## Number of obs: 66240, groups: Subject, 23
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.98254 0.10000 29.825 < 2e-16 ***
## Block.L 0.02700 0.05112 0.528 0.5974
## Block.Q -0.95787 0.05073 -18.883 < 2e-16 ***
## Block.C 0.26578 0.05313 5.002 5.66e-07 ***
## Block^4 -0.53388 0.05589 -9.553 < 2e-16 ***
## Block^5 0.12688 0.05700 2.226 0.0260 *
## Block^6 -0.13493 0.05901 -2.287 0.0222 *
## Block^7 0.11390 0.06000 1.898 0.0577 .
## Block^8 0.08343 0.05953 1.401 0.1611
## Block^9 -0.01478 0.06006 -0.246 0.8056
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Blck.L Blck.Q Blck.C Blck^4 Blck^5 Blck^6 Blck^7 Blck^8
## Block.L -0.004
## Block.Q -0.042 0.036
## Block.C 0.010 -0.298 0.026
## Block^4 -0.021 0.069 -0.218 0.005
## Block^5 0.005 -0.105 0.059 -0.122 -0.046
## Block^6 -0.005 0.038 -0.053 0.007 -0.088 -0.046
## Block^7 0.005 0.001 -0.001 -0.054 0.001 -0.027 -0.039
## Block^8 0.005 0.014 -0.029 -0.009 -0.017 0.001 0.006 -0.022
## Block^9 -0.001 0.011 0.014 -0.024 -0.017 0.015 0.017 -0.012 -0.038
# Obtain effects of the accuracy model
effects_acc <- allEffects(model_acc)
df_acc <- as.data.frame(effects_acc[[1]])
df_acc$l83 <- df_acc$fit - 1.3722 * df_acc$se
df_acc$u83 <- df_acc$fit + 1.3722 * df_acc$se
# Plot effects for Step Accuracy with error bars
ggplot(df_acc, aes(x = Block, y = fit)) +
geom_errorbar(aes(ymin = l83, ymax = u83), width = 0.2, size = 1, color = "aquamarine4") +
geom_point(color = "aquamarine3", size = 4) +
geom_vline(xintercept = 8.5, color = "red", size = 1) +
xlab("Block") + ylab("Step Accuracy (%)") +
ggtitle("B) Step.ACC ~ Block") +
theme_minimal() +
theme(text = element_text(size = base_size),
axis.title = element_text(size = base_size + 2),
plot.title = element_text(size = base_size + 4, face = "bold"),
axis.text = element_text(size = base_size - 2),
axis.text.y = element_text(face = "bold"),
axis.text.x = element_text(face = "bold"))
Anova(model_rt_lmm)
Anova(model_acc)
# Post hocs model RT
emm_blocks <- emmeans(model_rt_lmm, specs = ~ Block)
## Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
## To enable adjustments, add the argument 'pbkrtest.limit = 66240' (or larger)
## [or, globally, 'set emm_options(pbkrtest.limit = 66240)' or larger];
## but be warned that this may result in large computation time and memory use.
## Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
## To enable adjustments, add the argument 'lmerTest.limit = 66240' (or larger)
## [or, globally, 'set emm_options(lmerTest.limit = 66240)' or larger];
## but be warned that this may result in large computation time and memory use.
comparisons <- pairs(emm_blocks)
print(comparisons)
## contrast estimate SE df z.ratio p.value
## Block1 - Block2 327.84 11.8 Inf 27.763 <.0001
## Block1 - Block3 355.98 11.8 Inf 30.146 <.0001
## Block1 - Block4 380.00 11.8 Inf 32.181 <.0001
## Block1 - Block5 388.17 11.8 Inf 32.873 <.0001
## Block1 - Block6 388.86 11.8 Inf 32.931 <.0001
## Block1 - Block7 412.49 11.8 Inf 34.932 <.0001
## Block1 - Block8 431.55 11.8 Inf 36.546 <.0001
## Block1 - Block9 442.75 11.8 Inf 37.495 <.0001
## Block1 - Block10 366.19 11.8 Inf 31.011 <.0001
## Block2 - Block3 28.14 11.8 Inf 2.383 0.3360
## Block2 - Block4 52.16 11.8 Inf 4.417 0.0004
## Block2 - Block5 60.33 11.8 Inf 5.109 <.0001
## Block2 - Block6 61.02 11.8 Inf 5.167 <.0001
## Block2 - Block7 84.65 11.8 Inf 7.169 <.0001
## Block2 - Block8 103.71 11.8 Inf 8.783 <.0001
## Block2 - Block9 114.91 11.8 Inf 9.731 <.0001
## Block2 - Block10 38.35 11.8 Inf 3.248 0.0386
## Block3 - Block4 24.02 11.8 Inf 2.034 0.5745
## Block3 - Block5 32.19 11.8 Inf 2.726 0.1631
## Block3 - Block6 32.88 11.8 Inf 2.785 0.1416
## Block3 - Block7 56.51 11.8 Inf 4.786 0.0001
## Block3 - Block8 75.57 11.8 Inf 6.400 <.0001
## Block3 - Block9 86.77 11.8 Inf 7.348 <.0001
## Block3 - Block10 10.21 11.8 Inf 0.865 0.9974
## Block4 - Block5 8.17 11.8 Inf 0.692 0.9996
## Block4 - Block6 8.86 11.8 Inf 0.750 0.9992
## Block4 - Block7 32.49 11.8 Inf 2.752 0.1534
## Block4 - Block8 51.55 11.8 Inf 4.366 0.0005
## Block4 - Block9 62.75 11.8 Inf 5.314 <.0001
## Block4 - Block10 -13.81 11.8 Inf -1.169 0.9769
## Block5 - Block6 0.69 11.8 Inf 0.058 1.0000
## Block5 - Block7 24.32 11.8 Inf 2.060 0.5563
## Block5 - Block8 43.38 11.8 Inf 3.674 0.0090
## Block5 - Block9 54.58 11.8 Inf 4.622 0.0002
## Block5 - Block10 -21.98 11.8 Inf -1.861 0.6955
## Block6 - Block7 23.63 11.8 Inf 2.001 0.5981
## Block6 - Block8 42.69 11.8 Inf 3.615 0.0112
## Block6 - Block9 53.89 11.8 Inf 4.564 0.0002
## Block6 - Block10 -22.67 11.8 Inf -1.920 0.6556
## Block7 - Block8 19.06 11.8 Inf 1.614 0.8417
## Block7 - Block9 30.26 11.8 Inf 2.562 0.2357
## Block7 - Block10 -46.30 11.8 Inf -3.921 0.0035
## Block8 - Block9 11.20 11.8 Inf 0.948 0.9948
## Block8 - Block10 -65.36 11.8 Inf -5.535 <.0001
## Block9 - Block10 -76.56 11.8 Inf -6.483 <.0001
##
## Degrees-of-freedom method: asymptotic
## P value adjustment: tukey method for comparing a family of 10 estimates
# Post hocs model acc
emm_blocks <- emmeans(model_acc, specs = ~ Block)
comparisons <- pairs(emm_blocks)
print(comparisons)
## contrast estimate SE df z.ratio p.value
## Block1 - Block2 -1.08505 0.0722 Inf -15.034 <.0001
## Block1 - Block3 -1.12648 0.0731 Inf -15.403 <.0001
## Block1 - Block4 -1.06307 0.0717 Inf -14.832 <.0001
## Block1 - Block5 -1.13033 0.0732 Inf -15.434 <.0001
## Block1 - Block6 -1.02738 0.0709 Inf -14.486 <.0001
## Block1 - Block7 -0.94929 0.0693 Inf -13.705 <.0001
## Block1 - Block8 -1.05945 0.0716 Inf -14.798 <.0001
## Block1 - Block9 -0.83683 0.0670 Inf -12.491 <.0001
## Block1 - Block10 -0.33417 0.0590 Inf -5.665 <.0001
## Block2 - Block3 -0.04143 0.0866 Inf -0.478 1.0000
## Block2 - Block4 0.02197 0.0854 Inf 0.257 1.0000
## Block2 - Block5 -0.04528 0.0867 Inf -0.522 1.0000
## Block2 - Block6 0.05767 0.0848 Inf 0.680 0.9996
## Block2 - Block7 0.13576 0.0834 Inf 1.627 0.8352
## Block2 - Block8 0.02559 0.0854 Inf 0.300 1.0000
## Block2 - Block9 0.24822 0.0816 Inf 3.043 0.0711
## Block2 - Block10 0.75088 0.0751 Inf 9.992 <.0001
## Block3 - Block4 0.06341 0.0862 Inf 0.735 0.9993
## Block3 - Block5 -0.00385 0.0875 Inf -0.044 1.0000
## Block3 - Block6 0.09910 0.0856 Inf 1.157 0.9785
## Block3 - Block7 0.17720 0.0842 Inf 2.104 0.5247
## Block3 - Block8 0.06703 0.0862 Inf 0.778 0.9989
## Block3 - Block9 0.28965 0.0824 Inf 3.515 0.0160
## Block3 - Block10 0.79232 0.0761 Inf 10.418 <.0001
## Block4 - Block5 -0.06726 0.0863 Inf -0.779 0.9989
## Block4 - Block6 0.03569 0.0844 Inf 0.423 1.0000
## Block4 - Block7 0.11379 0.0830 Inf 1.371 0.9361
## Block4 - Block8 0.00362 0.0849 Inf 0.043 1.0000
## Block4 - Block9 0.22624 0.0811 Inf 2.789 0.1400
## Block4 - Block10 0.72891 0.0747 Inf 9.762 <.0001
## Block5 - Block6 0.10295 0.0857 Inf 1.201 0.9724
## Block5 - Block7 0.18104 0.0843 Inf 2.147 0.4941
## Block5 - Block8 0.07088 0.0862 Inf 0.822 0.9983
## Block5 - Block9 0.29350 0.0825 Inf 3.558 0.0137
## Block5 - Block10 0.79616 0.0762 Inf 10.452 <.0001
## Block6 - Block7 0.07809 0.0823 Inf 0.948 0.9948
## Block6 - Block8 -0.03207 0.0843 Inf -0.380 1.0000
## Block6 - Block9 0.19055 0.0805 Inf 2.368 0.3451
## Block6 - Block10 0.69321 0.0740 Inf 9.371 <.0001
## Block7 - Block8 -0.11017 0.0829 Inf -1.329 0.9472
## Block7 - Block9 0.11246 0.0790 Inf 1.423 0.9203
## Block7 - Block10 0.61512 0.0724 Inf 8.499 <.0001
## Block8 - Block9 0.22262 0.0811 Inf 2.746 0.1554
## Block8 - Block10 0.72529 0.0746 Inf 9.724 <.0001
## Block9 - Block10 0.50266 0.0702 Inf 7.161 <.0001
##
## Results are given on the log odds ratio (not the response) scale.
## P value adjustment: tukey method for comparing a family of 10 estimates