## Load execution‑phase data
setwd("C:/Users/Marcel/Desktop/Scripts")
df <- read_excel("hurst_df_xyz_v2.xlsx")
df$Block <- factor(df$Block, ordered = FALSE, levels = c('1','2','3','4','5','6','7','8','9','10'))
df$Trial <- as.numeric(df$Trial)
df$Subject <- factor(df$Subject)
df$trial_acc <- factor(df$trial_acc)
df_acc <- df %>%
filter(trial_acc == "1")
str(df)
## tibble [10,560 × 19] (S3: tbl_df/tbl/data.frame)
## $ Subject : Factor w/ 22 levels "11","12","13",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Block : Factor w/ 10 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Trial : num [1:10560] 1 2 3 4 5 6 7 8 9 10 ...
## $ Sequence : chr [1:10560] "1" "2" "1" "1" ...
## $ H_prep_X : num [1:10560] 0.867 0.996 0.894 0.719 0.843 ...
## $ H_exec_X : num [1:10560] 0.993 0.979 0.992 0.992 0.995 ...
## $ H_prep_Y : num [1:10560] 0.698 0.994 0.992 0.925 0.902 ...
## $ H_exec_Y : num [1:10560] 0.996 0.993 0.993 0.996 0.994 ...
## $ H_prep_Z : num [1:10560] 0.908 0.995 0.953 0.857 0.88 ...
## $ H_exec_Z : num [1:10560] 0.996 0.995 0.993 0.995 0.996 ...
## $ trial_mean_rt : num [1:10560] 1466 930 981 1210 1494 ...
## $ trial_total_rt: num [1:10560] 8796 5578 5884 7263 8965 ...
## $ trial_acc : Factor w/ 2 levels "0","1": 1 2 1 1 1 2 2 2 2 2 ...
## $ Logtrialmeanrt: num [1:10560] 7.29 6.83 6.89 7.1 7.31 ...
## $ Logtrialsumrt : num [1:10560] 9.08 8.63 8.68 8.89 9.1 ...
## $ Phase : chr [1:10560] "Learning" "Learning" "Learning" "Learning" ...
## $ expertise : num [1:10560] 0.834 0.834 0.834 0.834 0.834 ...
## $ H_prep_vm : num [1:10560] 0.817 0.995 0.987 0.84 0.797 ...
## $ H_exec_vm : num [1:10560] 0.994 0.993 0.991 0.996 0.993 ...
df$Block_Num <- as.numeric(as.character(df$Block))
# Subset for Learning and Test Phases
df_learning <- df %>% filter(Phase == "Learning")
df_test <- df %>% filter(Phase == "Test")
# Subset for Learning and Test Phases with only accurate Trials
df_learning_acc <- subset(df_acc, Block %in% c("1","2","3","4","5","6","7","8"))
df_test_acc <- subset(df_acc, Block %in% c("9", "10"))
df_learning_acc$Block <- factor(df_learning_acc$Block, levels = c('1','2','3','4','5','6','7','8'))
df_test_acc$Block <- factor(df_test_acc$Block, levels = c('9', '10'))
df_learning_acc$Subject <- factor(df_learning_acc$Subject)
df_test_acc$Subject <- factor(df_test_acc$Subject)
str(df_learning_acc)
## tibble [7,606 × 19] (S3: tbl_df/tbl/data.frame)
## $ Subject : Factor w/ 22 levels "11","12","13",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Block : Factor w/ 8 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Trial : num [1:7606] 2 6 7 8 9 10 12 13 14 15 ...
## $ Sequence : chr [1:7606] "2" "2" "2" "1" ...
## $ H_prep_X : num [1:7606] 0.996 0.951 0.989 0.988 0.954 ...
## $ H_exec_X : num [1:7606] 0.979 0.994 0.989 0.994 0.994 ...
## $ H_prep_Y : num [1:7606] 0.994 0.743 0.992 0.976 0.982 ...
## $ H_exec_Y : num [1:7606] 0.993 0.996 0.992 0.991 0.995 ...
## $ H_prep_Z : num [1:7606] 0.995 0.956 0.954 0.993 0.989 ...
## $ H_exec_Z : num [1:7606] 0.995 0.994 0.992 0.992 0.995 ...
## $ trial_mean_rt : num [1:7606] 930 1237 853 1031 837 ...
## $ trial_total_rt: num [1:7606] 5578 7420 5117 6188 5021 ...
## $ trial_acc : Factor w/ 2 levels "0","1": 2 2 2 2 2 2 2 2 2 2 ...
## $ Logtrialmeanrt: num [1:7606] 6.83 7.12 6.75 6.94 6.73 ...
## $ Logtrialsumrt : num [1:7606] 8.63 8.91 8.54 8.73 8.52 ...
## $ Phase : chr [1:7606] "Learning" "Learning" "Learning" "Learning" ...
## $ expertise : num [1:7606] 0.834 0.834 0.834 0.834 0.834 ...
## $ H_prep_vm : num [1:7606] 0.995 0.846 0.985 0.991 0.971 ...
## $ H_exec_vm : num [1:7606] 0.993 0.995 0.995 0.993 0.991 ...
str(df_test_acc)
## tibble [1,833 × 19] (S3: tbl_df/tbl/data.frame)
## $ Subject : Factor w/ 22 levels "11","12","13",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Block : Factor w/ 2 levels "9","10": 1 1 1 1 1 1 1 1 1 1 ...
## $ Trial : num [1:1833] 2 3 4 5 6 9 10 11 12 13 ...
## $ Sequence : chr [1:1833] "2" "1" "2" "2" ...
## $ H_prep_X : num [1:1833] 0.989 0.991 0.772 0.96 0.994 ...
## $ H_exec_X : num [1:1833] 0.974 0.938 0.976 0.977 0.973 ...
## $ H_prep_Y : num [1:1833] 0.985 0.99 0.808 0.938 0.994 ...
## $ H_exec_Y : num [1:1833] 0.908 0.937 0.973 0.973 0.982 ...
## $ H_prep_Z : num [1:1833] 0.976 0.991 0.822 0.913 0.987 ...
## $ H_exec_Z : num [1:1833] 0.994 0.978 0.991 0.991 0.994 ...
## $ trial_mean_rt : num [1:1833] 282 217 240 192 243 ...
## $ trial_total_rt: num [1:1833] 1689 1303 1437 1151 1458 ...
## $ trial_acc : Factor w/ 2 levels "0","1": 2 2 2 2 2 2 2 2 2 2 ...
## $ Logtrialmeanrt: num [1:1833] 5.64 5.38 5.48 5.26 5.49 ...
## $ Logtrialsumrt : num [1:1833] 7.43 7.17 7.27 7.05 7.28 ...
## $ Phase : chr [1:1833] "Test" "Test" "Test" "Test" ...
## $ expertise : num [1:1833] 0.834 0.834 0.834 0.834 0.834 ...
## $ H_prep_vm : num [1:1833] 0.979 0.994 0.715 0.956 0.993 ...
## $ H_exec_vm : num [1:1833] 0.972 0.904 0.975 0.969 0.981 ...
## ===============================================================
## 1) Learning phase model (Blocks 1–8) - Accuracy as DV
## ===============================================================
mod_acc_learn <- glmmTMB(
trial_acc ~ Block_Num * Logtrialmeanrt +
(1 | Subject),
family = binomial(),
data = df_learning
)
summary(mod_acc_learn)
## Family: binomial ( logit )
## Formula: trial_acc ~ Block_Num * Logtrialmeanrt + (1 | Subject)
## Data: df_learning
##
## AIC BIC logLik -2*log(L) df.resid
## 3718.6 3753.8 -1854.3 3708.6 8443
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Subject (Intercept) 2.246 1.499
## Number of obs: 8448, groups: Subject, 22
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 26.95565 1.26258 21.350 <2e-16 ***
## Block_Num 0.20653 0.24053 0.859 0.3905
## Logtrialmeanrt -3.84274 0.18752 -20.492 <2e-16 ***
## Block_Num:Logtrialmeanrt -0.06975 0.03844 -1.815 0.0696 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(mod_acc_learn)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: trial_acc
## Chisq Df Pr(>Chisq)
## Block_Num 114.2377 1 <2e-16 ***
## Logtrialmeanrt 938.0129 1 <2e-16 ***
## Block_Num:Logtrialmeanrt 3.2925 1 0.0696 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
p_acc_learn_2d <- plot_model(
mod_acc_learn,
type = "pred",
terms = c("Block_Num", "Logtrialmeanrt"),
transform = "response"
) +
labs(
x = "Block (numeric, 1–8)",
y = "Predicted accuracy",
title = "Predicted accuracy over practice (Learning phase)",
subtitle = "Lines = different log RT levels"
)
## You are calculating adjusted predictions on the population-level (i.e.
## `type = "fixed"`) for a *generalized* linear mixed model.
## This may produce biased estimates due to Jensen's inequality. Consider
## setting `bias_correction = TRUE` to correct for this bias.
## See also the documentation of the `bias_correction` argument.
print(p_acc_learn_2d)
## ===============================================================
## 2) Test phase model (Blocks 9–10) - Accuracy as DV
## ===============================================================
mod_acc_test <- glmmTMB(
trial_acc ~ Block_Num * Logtrialmeanrt +
(1 | Subject),
family = binomial(),
data = df_test
)
summary(mod_acc_test)
## Family: binomial ( logit )
## Formula: trial_acc ~ Block_Num * Logtrialmeanrt + (1 | Subject)
## Data: df_test
##
## AIC BIC logLik -2*log(L) df.resid
## 1177.6 1205.8 -583.8 1167.6 2107
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Subject (Intercept) 2.737 1.654
## Number of obs: 2112, groups: Subject, 22
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 43.0049 22.9818 1.871 0.0613 .
## Block_Num -1.3669 2.3965 -0.570 0.5684
## Logtrialmeanrt -7.2622 3.7688 -1.927 0.0540 .
## Block_Num:Logtrialmeanrt 0.2739 0.3924 0.698 0.4852
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(mod_acc_test)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: trial_acc
## Chisq Df Pr(>Chisq)
## Block_Num 3.1143 1 0.07761 .
## Logtrialmeanrt 271.4449 1 < 2e-16 ***
## Block_Num:Logtrialmeanrt 0.4872 1 0.48516
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
p_acc_test_2d <- plot_model(
mod_acc_test,
type = "pred",
terms = c("Block_Num", "Logtrialmeanrt"),
transform = "response"
) +
labs(
x = "Block (numeric, 9–10)",
y = "Predicted accuracy",
title = "Predicted accuracy over practice (Test phase)",
subtitle = "Lines = different log RT levels"
)
print(p_acc_test_2d)
## ===============================================================
## 1) Learning phase model (Blocks 1–8) - Logtrialmeanrt as DV
## ===============================================================
model_learning <- glmmTMB(
Logtrialmeanrt ~ H_prep_vm * H_exec_vm * Block +
(1 | Subject),
data = df_learning_acc,
family = gaussian()
)
summary(model_learning)
## Family: gaussian ( identity )
## Formula:
## Logtrialmeanrt ~ H_prep_vm * H_exec_vm * Block + (1 | Subject)
## Data: df_learning_acc
##
## AIC BIC logLik -2*log(L) df.resid
## 1687.9 1923.7 -809.9 1619.9 7572
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Subject (Intercept) 0.11729 0.3425
## Residual 0.07113 0.2667
## Number of obs: 7606, groups: Subject, 22
##
## Dispersion estimate for gaussian family (sigma^2): 0.0711
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 25.525 4.611 5.536 3.10e-08 ***
## H_prep_vm -25.449 5.200 -4.894 9.86e-07 ***
## H_exec_vm -19.506 4.650 -4.195 2.73e-05 ***
## Block2 -19.482 5.387 -3.617 0.000298 ***
## Block3 -17.232 5.593 -3.081 0.002061 **
## Block4 -16.461 5.251 -3.135 0.001718 **
## Block5 -16.275 5.389 -3.020 0.002526 **
## Block6 -22.082 5.537 -3.988 6.67e-05 ***
## Block7 -22.889 5.314 -4.307 1.65e-05 ***
## Block8 -15.663 5.147 -3.043 0.002340 **
## H_prep_vm:H_exec_vm 25.902 5.245 4.938 7.88e-07 ***
## H_prep_vm:Block2 23.049 6.214 3.709 0.000208 ***
## H_prep_vm:Block3 20.660 6.316 3.271 0.001071 **
## H_prep_vm:Block4 19.650 5.992 3.279 0.001041 **
## H_prep_vm:Block5 20.796 6.096 3.411 0.000647 ***
## H_prep_vm:Block6 26.685 6.281 4.248 2.15e-05 ***
## H_prep_vm:Block7 27.454 5.988 4.584 4.55e-06 ***
## H_prep_vm:Block8 18.049 5.811 3.106 0.001898 **
## H_exec_vm:Block2 19.529 5.435 3.593 0.000327 ***
## H_exec_vm:Block3 17.342 5.644 3.072 0.002123 **
## H_exec_vm:Block4 16.386 5.300 3.091 0.001992 **
## H_exec_vm:Block5 16.106 5.440 2.960 0.003072 **
## H_exec_vm:Block6 21.906 5.591 3.918 8.93e-05 ***
## H_exec_vm:Block7 22.704 5.366 4.231 2.33e-05 ***
## H_exec_vm:Block8 15.229 5.196 2.931 0.003379 **
## H_prep_vm:H_exec_vm:Block2 -23.519 6.272 -3.750 0.000177 ***
## H_prep_vm:H_exec_vm:Block3 -21.286 6.376 -3.338 0.000843 ***
## H_prep_vm:H_exec_vm:Block4 -20.131 6.051 -3.327 0.000878 ***
## H_prep_vm:H_exec_vm:Block5 -21.238 6.156 -3.450 0.000561 ***
## H_prep_vm:H_exec_vm:Block6 -27.131 6.344 -4.277 1.90e-05 ***
## H_prep_vm:H_exec_vm:Block7 -27.959 6.049 -4.622 3.79e-06 ***
## H_prep_vm:H_exec_vm:Block8 -18.269 5.868 -3.113 0.001851 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(model_learning)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Logtrialmeanrt
## Chisq Df Pr(>Chisq)
## H_prep_vm 0.3157 1 0.5741815
## H_exec_vm 59.5933 1 1.166e-14 ***
## Block 3220.1052 7 < 2.2e-16 ***
## H_prep_vm:H_exec_vm 14.6604 1 0.0001287 ***
## H_prep_vm:Block 38.1020 7 2.898e-06 ***
## H_exec_vm:Block 17.1039 7 0.0167384 *
## H_prep_vm:H_exec_vm:Block 25.9254 7 0.0005194 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot_model(
model_learning, type = "pred",
terms = c("H_exec_vm", "H_prep_vm", "Block")
)
# Block Pairwise comparisons
emm_learning <- emmeans(model_learning, ~ Block)
## NOTE: Results may be misleading due to involvement in interactions
print(emm_learning)
## Block emmean SE df lower.CL upper.CL
## 1 6.37 0.0737 7572 6.22 6.51
## 2 6.02 0.0735 7572 5.88 6.17
## 3 5.96 0.0735 7572 5.82 6.10
## 4 5.90 0.0735 7572 5.75 6.04
## 5 5.85 0.0735 7572 5.71 6.00
## 6 5.83 0.0735 7572 5.68 5.97
## 7 5.77 0.0735 7572 5.63 5.91
## 8 5.75 0.0735 7572 5.61 5.89
##
## Confidence level used: 0.95
# All pairwise Block comparisons (1–8)
contr_learning <- pairs(emm_learning)
print(contr_learning)
## contrast estimate SE df t.ratio p.value
## Block1 - Block2 0.3446 0.0135 7572 25.621 <.0001
## Block1 - Block3 0.4085 0.0134 7572 30.407 <.0001
## Block1 - Block4 0.4723 0.0134 7572 35.185 <.0001
## Block1 - Block5 0.5167 0.0134 7572 38.475 <.0001
## Block1 - Block6 0.5393 0.0135 7572 40.039 <.0001
## Block1 - Block7 0.5982 0.0136 7572 44.115 <.0001
## Block1 - Block8 0.6178 0.0136 7572 45.385 <.0001
## Block2 - Block3 0.0639 0.0123 7572 5.211 <.0001
## Block2 - Block4 0.1277 0.0123 7572 10.423 <.0001
## Block2 - Block5 0.1721 0.0123 7572 14.032 <.0001
## Block2 - Block6 0.1947 0.0123 7572 15.840 <.0001
## Block2 - Block7 0.2535 0.0124 7572 20.486 <.0001
## Block2 - Block8 0.2732 0.0124 7572 22.004 <.0001
## Block3 - Block4 0.0638 0.0121 7572 5.257 <.0001
## Block3 - Block5 0.1082 0.0122 7572 8.901 <.0001
## Block3 - Block6 0.1308 0.0122 7572 10.745 <.0001
## Block3 - Block7 0.1896 0.0122 7572 15.483 <.0001
## Block3 - Block8 0.2093 0.0123 7572 17.049 <.0001
## Block4 - Block5 0.0444 0.0121 7572 3.654 0.0063
## Block4 - Block6 0.0670 0.0122 7572 5.506 <.0001
## Block4 - Block7 0.1258 0.0122 7572 10.282 <.0001
## Block4 - Block8 0.1455 0.0123 7572 11.865 <.0001
## Block5 - Block6 0.0226 0.0122 7572 1.854 0.5832
## Block5 - Block7 0.0814 0.0123 7572 6.644 <.0001
## Block5 - Block8 0.1011 0.0123 7572 8.231 <.0001
## Block6 - Block7 0.0589 0.0123 7572 4.800 <.0001
## Block6 - Block8 0.0785 0.0123 7572 6.393 <.0001
## Block7 - Block8 0.0197 0.0124 7572 1.593 0.7548
##
## P value adjustment: tukey method for comparing a family of 8 estimates
print(confint(contr_learning))
## contrast estimate SE df lower.CL upper.CL
## Block1 - Block2 0.3446 0.0135 7572 0.30385 0.3854
## Block1 - Block3 0.4085 0.0134 7572 0.36779 0.4492
## Block1 - Block4 0.4723 0.0134 7572 0.43164 0.5130
## Block1 - Block5 0.5167 0.0134 7572 0.47600 0.5574
## Block1 - Block6 0.5393 0.0135 7572 0.49846 0.5801
## Block1 - Block7 0.5982 0.0136 7572 0.55705 0.6393
## Block1 - Block8 0.6178 0.0136 7572 0.57657 0.6591
## Block2 - Block3 0.0639 0.0123 7572 0.02672 0.1011
## Block2 - Block4 0.1277 0.0123 7572 0.09056 0.1649
## Block2 - Block5 0.1721 0.0123 7572 0.13491 0.2093
## Block2 - Block6 0.1947 0.0123 7572 0.15741 0.2319
## Block2 - Block7 0.2535 0.0124 7572 0.21601 0.2911
## Block2 - Block8 0.2732 0.0124 7572 0.23557 0.3109
## Block3 - Block4 0.0638 0.0121 7572 0.02701 0.1006
## Block3 - Block5 0.1082 0.0122 7572 0.07135 0.1451
## Block3 - Block6 0.1308 0.0122 7572 0.09388 0.1677
## Block3 - Block7 0.1896 0.0122 7572 0.15251 0.2268
## Block3 - Block8 0.2093 0.0123 7572 0.17210 0.2465
## Block4 - Block5 0.0444 0.0121 7572 0.00756 0.0812
## Block4 - Block6 0.0670 0.0122 7572 0.03009 0.1038
## Block4 - Block7 0.1258 0.0122 7572 0.08872 0.1629
## Block4 - Block8 0.1455 0.0123 7572 0.10832 0.1827
## Block5 - Block6 0.0226 0.0122 7572 -0.01435 0.0595
## Block5 - Block7 0.0814 0.0123 7572 0.04428 0.1186
## Block5 - Block8 0.1011 0.0123 7572 0.06387 0.1384
## Block6 - Block7 0.0589 0.0123 7572 0.02169 0.0960
## Block6 - Block8 0.0785 0.0123 7572 0.04129 0.1158
## Block7 - Block8 0.0197 0.0124 7572 -0.01777 0.0571
##
## Confidence level used: 0.95
## Conf-level adjustment: tukey method for comparing a family of 8 estimates
## ===============================================================
## 2) Test phase model (Blocks 9–10) - Logtrialmeanrt as DV
## ===============================================================
model_test <- glmmTMB(
Logtrialmeanrt ~ H_prep_vm * H_exec_vm * Block +
(1 | Subject),
family = gaussian(),
data = df_test_acc
)
summary(model_test)
## Family: gaussian ( identity )
## Formula:
## Logtrialmeanrt ~ H_prep_vm * H_exec_vm * Block + (1 | Subject)
## Data: df_test_acc
##
## AIC BIC logLik -2*log(L) df.resid
## 224.9 280.1 -102.5 204.9 1823
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Subject (Intercept) 0.13018 0.3608
## Residual 0.06153 0.2481
## Number of obs: 1833, groups: Subject, 22
##
## Dispersion estimate for gaussian family (sigma^2): 0.0615
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.63935 3.33337 1.692 0.0907 .
## H_prep_vm -2.65272 3.57201 -0.743 0.4577
## H_exec_vm -0.09203 3.37772 -0.027 0.9783
## Block10 -6.90770 4.27426 -1.616 0.1061
## H_prep_vm:H_exec_vm 2.87813 3.62177 0.795 0.4268
## H_prep_vm:Block10 6.27672 4.66034 1.347 0.1780
## H_exec_vm:Block10 7.53499 4.32775 1.741 0.0817 .
## H_prep_vm:H_exec_vm:Block10 -6.74790 4.72071 -1.429 0.1529
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(model_test)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Logtrialmeanrt
## Chisq Df Pr(>Chisq)
## H_prep_vm 0.0007 1 0.979249
## H_exec_vm 74.6506 1 < 2.2e-16 ***
## Block 289.5205 1 < 2.2e-16 ***
## H_prep_vm:H_exec_vm 0.1754 1 0.675315
## H_prep_vm:Block 22.8309 1 1.769e-06 ***
## H_exec_vm:Block 8.7220 1 0.003144 **
## H_prep_vm:H_exec_vm:Block 2.0433 1 0.152882
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
p_test <- plot_model(
model_test,
type = "pred",
terms = c("H_exec_vm", "H_prep_vm", "Block")
)
p_test + labs(
title = "Predicted log RT surfaces",
subtitle = "Separate surfaces for Block 9 (familiar) and Block 10 (unfamiliar)",
x = "H_exec_vm",
y = "H_prep_vm"
)
# Block Pairwise comparisons
emm_test <- emmeans(model_test, ~ Block)
## NOTE: Results may be misleading due to involvement in interactions
print(emm_test)
## Block emmean SE df lower.CL upper.CL
## 9 5.70 0.0774 1823 5.54 5.85
## 10 5.89 0.0774 1823 5.74 6.04
##
## Confidence level used: 0.95
# All pairwise Block comparisons (1–8)
contr_test <- pairs(emm_test)
print(contr_test)
## contrast estimate SE df t.ratio p.value
## Block9 - Block10 -0.194 0.0119 1823 -16.308 <.0001
print(confint(contr_test))
## contrast estimate SE df lower.CL upper.CL
## Block9 - Block10 -0.194 0.0119 1823 -0.218 -0.171
##
## Confidence level used: 0.95
## ===============================================================
## STITCHED SURFACES: Blocks 1–8 from model_learning,
## Blocks 9–10 from model_test
## ===============================================================
## ---------------------------------------------------------------
## Hyperparameters
## ---------------------------------------------------------------
q_lower <- 0.025
q_upper <- 0.975
q_lower_pct <- q_lower * 100
q_upper_pct <- q_upper * 100
viridis_option <- "viridis"
viridis_begin <- 0
viridis_end <- 1
fast_color <- "light"
viridis_direction <- if (fast_color == "dark") 1L else -1L
# Plot theme
common_theme <- theme_classic(base_size = 14) +
theme(
axis.title.x = element_text(size = 14),
axis.title.y = element_text(size = 14),
axis.text.x = element_text(size = 9),
axis.text.y = element_text(size = 9),
strip.text = element_text(size = 12, face = "bold"),
strip.background = element_rect(fill = "grey90", color = NA),
legend.title = element_text(size = 12),
legend.text = element_text(size = 10)
)
## ---------------------------------------------------------------
## Block sets and labels
## ---------------------------------------------------------------
blocks_learning <- levels(df_learning_acc$Block)
blocks_test <- levels(df_test_acc$Block)
blocks_all <- c(blocks_learning, blocks_test)
block_labs_all <- setNames(
paste0("Block ", blocks_all),
blocks_all
)
# block_labs_all["10"] <- "Block 10 (unfamiliar)"
## ===============================================================
## A. WITHIN-BLOCK SURFACES (per-block H ranges)
## ===============================================================
## ---------------------------
## Learning within-block grids
## ---------------------------
grid_list_block_learning <- lapply(blocks_learning, function(b) {
df_b <- df_learning_acc %>% filter(Block == b)
H_exec_seq_b <- seq(
quantile(df_b$H_exec_vm, q_lower, na.rm = TRUE),
quantile(df_b$H_exec_vm, q_upper, na.rm = TRUE),
length.out = 100
)
H_prep_seq_b <- seq(
quantile(df_b$H_prep_vm, q_lower, na.rm = TRUE),
quantile(df_b$H_prep_vm, q_upper, na.rm = TRUE),
length.out = 100
)
expand.grid(
H_exec_vm = H_exec_seq_b,
H_prep_vm = H_prep_seq_b,
Block = b,
KEEP.OUT.ATTRS = FALSE,
stringsAsFactors = FALSE
)
})
pred_grid_block_learning <- bind_rows(grid_list_block_learning)
pred_grid_block_learning$Block <- factor(
pred_grid_block_learning$Block,
levels = blocks_all
)
pred_grid_block_learning <- as.data.frame(pred_grid_block_learning)
pred_grid_block_learning$pred_logrt <- predict(
model_learning,
newdata = pred_grid_block_learning,
type = "response",
re.form = NA,
allow.new.levels = TRUE
)
## ---------------------------
## Test within-block grids
## ---------------------------
grid_list_block_test <- lapply(blocks_test, function(b) {
df_b <- df_test_acc %>% filter(Block == b)
H_exec_seq_b <- seq(
quantile(df_b$H_exec_vm, q_lower, na.rm = TRUE),
quantile(df_b$H_exec_vm, q_upper, na.rm = TRUE),
length.out = 100
)
H_prep_seq_b <- seq(
quantile(df_b$H_prep_vm, q_lower, na.rm = TRUE),
quantile(df_b$H_prep_vm, q_upper, na.rm = TRUE),
length.out = 100
)
expand.grid(
H_exec_vm = H_exec_seq_b,
H_prep_vm = H_prep_seq_b,
Block = b,
KEEP.OUT.ATTRS = FALSE,
stringsAsFactors = FALSE
)
})
pred_grid_block_test <- bind_rows(grid_list_block_test)
pred_grid_block_test$Block <- factor(
pred_grid_block_test$Block,
levels = blocks_all
)
pred_grid_block_test <- as.data.frame(pred_grid_block_test)
pred_grid_block_test$pred_logrt <- predict(
model_test,
newdata = pred_grid_block_test,
type = "response",
re.form = NA,
allow.new.levels = TRUE
)
## ---------------------------
## Combine within-block grids
## ---------------------------
pred_grid_block_all <- bind_rows(
pred_grid_block_learning,
pred_grid_block_test
)
## ===============================================================
## B. COMMON COLOR RANGE AND THEME FOR MODEL SURFACES
## ===============================================================
rt_range_all <- range(pred_grid_block_all$pred_logrt, na.rm = TRUE)
## ===============================================================
## C. STITCHED MODEL HEATMAPS (WITHIN-BLOCK RANGES)
## ===============================================================
p_block_all <- ggplot() +
geom_tile(
data = pred_grid_block_all,
aes(
x = H_prep_vm,
y = H_exec_vm,
fill = pred_logrt
)
) +
facet_wrap(
~ Block,
ncol = 5,
labeller = labeller(Block = block_labs_all)
) +
scale_fill_viridis_c(
option = viridis_option,
begin = viridis_begin,
end = viridis_end,
direction = viridis_direction,
limits = rt_range_all,
name = "Predicted\nlog RT"
) +
labs(
title = "Predicted RT surfaces",
subtitle = sprintf(
"H ranges = block-wise %.1f–%.1f%% quantiles",
q_lower_pct, q_upper_pct
),
x = "Preparation H",
y = "Execution H"
) +
common_theme
print(p_block_all)
## ===============================================================
## D. STITCHED RAW-DATA HEATMAPS
## ===============================================================
df_all <- bind_rows(
df_learning_acc %>% mutate(phase = "learning"),
df_test_acc %>% mutate(phase = "test")
)
df_all$Block <- factor(df_all$Block, levels = blocks_all)
rt_range_raw <- range(df_all$Logtrialmeanrt, na.rm = TRUE)
p_raw_all <- ggplot(
df_all,
aes(x = H_prep_vm, y = H_exec_vm)
) +
stat_summary_2d(
aes(z = Logtrialmeanrt),
bins = 30,
fun = mean
) +
facet_wrap(
~ Block,
ncol = 5,
labeller = labeller(Block = block_labs_all)
) +
scale_fill_viridis_c(
option = viridis_option,
begin = viridis_begin,
end = viridis_end,
direction = viridis_direction,
limits = rt_range_raw,
name = "Mean\nlog RT"
) +
labs(
title = "Raw-data RT heatmaps (Blocks 1–10)",
subtitle = "Mean log RT in H_prep × H_exec bins",
x = "Preparation H",
y = "Execution H"
) +
common_theme
print(p_raw_all)
summary(df_acc$H_exec_vm)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.7440 0.9875 0.9917 0.9851 0.9938 0.9972
summary(df_acc$H_prep_vm)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.3573 0.7234 0.8875 0.8419 0.9816 0.9973
## ===============================================================
## H prep / H exec over numeric Blocks for Learning and Test
## ===============================================================
df_learning_acc$Block_Num <- as.numeric(as.character(df_learning_acc$Block))
df_test_acc$Block_Num <- as.numeric(as.character(df_test_acc$Block))
# H_exec over blocks LEARNING
mod_Hexec_learning <- glmmTMB(
H_exec_vm ~ Block_Num +
(1 | Subject),
data = df_learning_acc,
family = gaussian()
)
# H_prep over blocks LEARNING
mod_Hprep_learning <- glmmTMB(
H_prep_vm ~ Block_Num +
(1 | Subject),
data = df_learning_acc,
family = gaussian()
)
plot_model(mod_Hexec_learning, type = "pred", terms = "Block_Num")
plot_model(mod_Hprep_learning, type = "pred", terms = "Block_Num")
summary(mod_Hexec_learning)
## Family: gaussian ( identity )
## Formula: H_exec_vm ~ Block_Num + (1 | Subject)
## Data: df_learning_acc
##
## AIC BIC logLik -2*log(L) df.resid
## -42757.7 -42730.0 21382.9 -42765.7 7602
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Subject (Intercept) 0.0002555 0.01598
## Residual 0.0002080 0.01442
## Number of obs: 7606, groups: Subject, 22
##
## Dispersion estimate for gaussian family (sigma^2): 0.000208
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 9.907e-01 3.428e-03 288.99 <2e-16 ***
## Block_Num -1.112e-03 7.316e-05 -15.19 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(mod_Hprep_learning)
## Family: gaussian ( identity )
## Formula: H_prep_vm ~ Block_Num + (1 | Subject)
## Data: df_learning_acc
##
## AIC BIC logLik -2*log(L) df.resid
## -9575.1 -9547.4 4791.5 -9583.1 7602
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Subject (Intercept) 0.005513 0.07425
## Residual 0.016382 0.12799
## Number of obs: 7606, groups: Subject, 22
##
## Dispersion estimate for gaussian family (sigma^2): 0.0164
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.8160476 0.0161705 50.47 <2e-16 ***
## Block_Num 0.0056147 0.0006492 8.65 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# -----------------------------------------------------------------------------------
# H_exec over blocks TEST
mod_Hexec_test <- glmmTMB(
H_exec_vm ~ Block_Num +
(1 | Subject),
data = df_test_acc,
family = gaussian()
)
# H_prep over blocks TEST
mod_Hprep_test <- glmmTMB(
H_prep_vm ~ Block_Num +
(1 | Subject),
data = df_test_acc,
family = gaussian()
)
plot_model(mod_Hexec_test, type = "pred", terms = "Block_Num")
plot_model(mod_Hprep_test, type = "pred", terms = "Block_Num")
summary(mod_Hexec_test)
## Family: gaussian ( identity )
## Formula: H_exec_vm ~ Block_Num + (1 | Subject)
## Data: df_test_acc
##
## AIC BIC logLik -2*log(L) df.resid
## -9957.1 -9935.0 4982.5 -9965.1 1829
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Subject (Intercept) 0.0003976 0.01994
## Residual 0.0002403 0.01550
## Number of obs: 1833, groups: Subject, 22
##
## Dispersion estimate for gaussian family (sigma^2): 0.00024
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.9733431 0.0081009 120.15 <2e-16 ***
## Block_Num 0.0008839 0.0007261 1.22 0.223
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(mod_Hprep_test)
## Family: gaussian ( identity )
## Formula: H_prep_vm ~ Block_Num + (1 | Subject)
## Data: df_test_acc
##
## AIC BIC logLik -2*log(L) df.resid
## -2503.9 -2481.8 1255.9 -2511.9 1829
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Subject (Intercept) 0.007384 0.08593
## Residual 0.014212 0.11921
## Number of obs: 1833, groups: Subject, 22
##
## Dispersion estimate for gaussian family (sigma^2): 0.0142
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.788087 0.056103 14.047 <2e-16 ***
## Block_Num 0.006705 0.005583 1.201 0.23
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## ---------- EXECUTION H: familiar (1–9) + unfamiliar (10) ----------
p_exec <- ggplot(df_acc, aes(x = Block, y = H_exec_vm)) +
geom_boxplot(fill = "aquamarine4", alpha = 0.35, outlier.shape = NA) +
geom_jitter(width = 0.15, alpha = 0.4, size = 1, colour = "aquamarine4") +
# vertical line between Block 9 and 10
geom_vline(xintercept = 9.5, linetype = "solid", colour = "aquamarine4", linewidth = 1) +
labs(x = "Block", y = "H (execution)",
title = "Execution H – familiar (1–9) vs unfamiliar (10)") +
theme_classic()
## ---------- PREPARATION H: familiar (1–9) + unfamiliar (10) ----------
p_prep <- ggplot(df_acc, aes(x = Block, y = H_prep_vm)) +
geom_boxplot(fill = "darkorchid", alpha = 0.35, outlier.shape = NA) +
geom_jitter(width = 0.15, alpha = 0.4, size = 1, colour = "darkorchid") +
geom_vline(xintercept = 9.5, linetype = "solid", colour = "darkorchid", linewidth = 1) +
labs(x = "Block", y = "H (preparation)",
title = "Preparation H – familiar (1–9) vs unfamiliar (10)") +
theme_classic()
p_exec
p_prep
## ---------- EXECUTION H: familiar (1–9) + unfamiliar (10) capped ----------
p_exec <- ggplot(df_acc, aes(x = Block, y = H_exec_vm)) +
geom_boxplot(fill = "aquamarine4", alpha = 0.35, outlier.shape = NA) +
geom_jitter(width = 0.15, alpha = 0.4, size = 1, colour = "aquamarine4") +
# vertical line between Block 9 and 10
geom_vline(xintercept = 9.5, linetype = "solid", colour = "aquamarine4", linewidth = 1) +
scale_y_continuous(limits = c(0.985, 1.0)) +
labs(x = "Block", y = "H (execution)",
title = "Execution H – familiar (1–9) vs unfamiliar (10)") +
theme_classic()
## ---------- PREPARATION H: familiar (1–9) + unfamiliar (10) capped----------
p_prep <- ggplot(df_acc, aes(x = Block, y = H_prep_vm)) +
geom_boxplot(fill = "darkorchid", alpha = 0.35, outlier.shape = NA) +
geom_jitter(width = 0.15, alpha = 0.4, size = 1, colour = "darkorchid") +
geom_vline(xintercept = 9.5, linetype = "solid", colour = "darkorchid", linewidth = 1) +
scale_y_continuous(limits = c(0.6, 1.0)) +
labs(x = "Block", y = "H (preparation)",
title = "Preparation H – familiar (1–9) vs unfamiliar (10)") +
theme_classic()
p_exec
## Warning: Removed 1800 rows containing non-finite outside the scale range
## (`stat_boxplot()`).
## Warning: Removed 1800 rows containing missing values or values outside the scale range
## (`geom_point()`).
p_prep
## Warning: Removed 780 rows containing non-finite outside the scale range
## (`stat_boxplot()`).
## Warning: Removed 780 rows containing missing values or values outside the scale range
## (`geom_point()`).
## ===============================================================
## Mosaic of prep/exec H quadrants by block
## - Area shows number of trials
## - Fill shows mean log RT
## ===============================================================
## --- Global bins across all 10 blocks ---------------------------
prep_breaks <- quantile(df_acc$H_prep_vm, probs = c(0, 0.5, 1), na.rm = TRUE)
exec_breaks <- quantile(df_acc$H_exec_vm, probs = c(0, 0.5, 1), na.rm = TRUE)
df_quad <- df_acc %>%
mutate(
prep_bin = cut(
H_prep_vm,
breaks = prep_breaks,
labels = c("🡻 prep H", " 🢁 prep H"),
include.lowest = TRUE
),
exec_bin = cut(
H_exec_vm,
breaks = exec_breaks,
labels = c("🢀 exec H", " 🢂 exec H"),
include.lowest = TRUE
)
) %>%
group_by(Block, prep_bin, exec_bin) %>%
summarise(
mean_rt = mean(Logtrialmeanrt, na.rm = TRUE),
n = n(),
.groups = "drop"
)
## --- Mosaic coordinates (area proportional to n) -----------------
prep_marg <- df_quad %>%
group_by(Block, prep_bin) %>%
summarise(n_prep = sum(n), .groups = "drop") %>%
group_by(Block) %>%
mutate(
total_n = sum(n_prep),
width = n_prep / total_n,
x_min = dplyr::lag(cumsum(width), default = 0),
x_max = cumsum(width)
) %>%
ungroup()
df_mosaic <- df_quad %>%
left_join(prep_marg, by = c("Block", "prep_bin")) %>%
group_by(Block, prep_bin) %>%
mutate(
height = n / n_prep,
y_min = dplyr::lag(cumsum(height), default = 0),
y_max = cumsum(height)
) %>%
ungroup()
## --- Axis label positions (tighter) ------------------------------------
prep_labels <- prep_marg %>%
mutate(
x = (x_min + x_max) / 2,
y = -0.06
)
exec_marg2 <- df_mosaic %>%
group_by(Block, exec_bin) %>%
summarise(
y_min = min(y_min),
y_max = max(y_max),
.groups = "drop"
) %>%
mutate(
y = (y_min + y_max) / 2,
x = -0.06
)
## --- Facet layout info -----------------------------------------------
n_cols <- 5
n_blocks <- length(levels(df_acc$Block))
n_rows <- ceiling(n_blocks / n_cols)
block_levels <- levels(df_acc$Block)
bottom_row_blocks <- block_levels[((n_rows - 1) * n_cols + 1):n_blocks]
left_col_blocks <- block_levels[seq(1, n_blocks, by = n_cols)]
## --- Plot --------------------------------------------------------------
gg_mosaic <- ggplot() +
geom_rect(
data = df_mosaic,
aes(
xmin = x_min,
xmax = x_max,
ymin = y_min,
ymax = y_max,
fill = mean_rt
),
color = "grey20"
) +
geom_text(
data = prep_labels %>% filter(Block %in% bottom_row_blocks),
aes(x = x, y = y, label = prep_bin),
inherit.aes = FALSE,
size = 3
) +
geom_text(
data = exec_marg2 %>% filter(Block %in% left_col_blocks),
aes(x = x, y = y, label = exec_bin),
inherit.aes = FALSE,
angle = 90,
hjust = 0.5,
size = 3
) +
facet_wrap(
~ Block,
nrow = 2,
labeller = labeller(Block = block_labs_all)
) +
scale_fill_viridis_c(
option = "viridis",
direction = -1,
name = "Mean\nlog RT"
) +
coord_equal(clip = "off") +
labs(
x = NULL,
y = NULL,
title = "Prep/exec H quadrants over blocks",
subtitle = "Bins defined by global median splits; quadrant area reflects trial count, color indicates mean log RT"
) +
theme_minimal() +
theme(
panel.grid = element_blank(),
axis.text = element_blank(),
axis.ticks = element_blank(),
strip.background = element_rect(fill = "grey90", color = NA),
strip.text = element_text(size = 12, face = "bold"),
plot.margin = margin(10, 10, 25, 25),
legend.position = "right",
legend.margin = margin(t = 2, r = 0, b = 0, l = 0)
)
gg_mosaic
df_acc$Block_Num <- as.numeric(as.character(df_acc$Block))
# Create a binary variable that is 1 if the Trial is in the High prep Low Exec Bin, 0 otherwise
df_acc <- df_acc %>%
mutate(
high_prep = H_prep_vm > median(H_prep_vm, na.rm = TRUE),
low_exec = H_exec_vm < median(H_exec_vm, na.rm = TRUE),
highprep_lowexec = as.integer(high_prep & low_exec)
)
# Modelling probability of observing a Trial within the High prep low exec bin over Blocks - with Block as NUMERIC
# positive significant slope indicates that over the blocks more trials populate that bin
mod_quad_num <- glmmTMB(
highprep_lowexec ~ Block_Num +
(1 | Subject),
family = binomial(),
data = df_acc
)
summary(mod_quad_num)
## Family: binomial ( logit )
## Formula: highprep_lowexec ~ Block_Num + (1 | Subject)
## Data: df_acc
##
## AIC BIC logLik -2*log(L) df.resid
## 8995.6 9017.0 -4494.8 8989.6 9436
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Subject (Intercept) 2.105 1.451
## Number of obs: 9439, groups: Subject, 22
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.12589 0.31698 -6.707 1.99e-11 ***
## Block_Num 0.12550 0.00933 13.452 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(mod_quad_num)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: highprep_lowexec
## Chisq Df Pr(>Chisq)
## Block_Num 180.95 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1