Lyapunov_Exponent_Analysis

library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.1     ✔ stringr   1.5.2
## ✔ ggplot2   4.0.0     ✔ tibble    3.3.0
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.1.0     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(lme4)

## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## 
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack

library(lmerTest)

## 
## Attaching package: 'lmerTest'
## 
## The following object is masked from 'package:lme4':
## 
##     lmer
## 
## The following object is masked from 'package:stats':
## 
##     step

library(ggplot2)
library(ggeffects)
library(dplyr)
library(patchwork)
library(emmeans)

## Welcome to emmeans.
## Caution: You lose important information if you filter this package's results.
## See '? untidy'

library(scales)

## 
## Attaching package: 'scales'
## 
## The following object is masked from 'package:purrr':
## 
##     discard
## 
## The following object is masked from 'package:readr':
## 
##     col_factor

library(cluster)
library(factoextra)

## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa

library(tidyr)
library(broom)
library(plotly)

## 
## Attaching package: 'plotly'
## 
## 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

# Configuration (must match Python script settings)
PARTICIPANT_IDS <- c(2, 3, 4, 5, 7, 8, 10, 11, 13, 14, 15, 16, 17, 18, 19, 20, 22, 23) # Placeholder for the participant IDs being analyzed
BASE_XSENS_DIR <- "/Users/vivi/UT/Thesis/Data/Xsens/Participants"
LYE_FAILURE_SENTINEL <- -9999.0 # Value used by the Python script for failed LLE calculations

load_lyap_data <- function(pid, period) {
  period_to_folder_map <- list(
    "prep" = "Movement Preparation Period",
    "exec" = "Movement Execution Period"
  )
  period_to_name_map <- list(
    "prep" = "Prep",
    "exec" = "Exec"
  )
  
  lyap_col_name <- paste0("Lyap_", period)
  
  folder_name <- period_to_folder_map[[period]]
  name_part <- period_to_name_map[[period]]
  
  file_name <- paste0("ID", pid, "_", name_part, "_LyE.csv")
  file_path <- file.path(BASE_XSENS_DIR, folder_name, file_name)
  
  if (!file.exists(file_path)) {
    stop(paste("Error: Data file not found:", file_path))
  }
  
  #message(paste("Loading data from:", file_path))
  
  df <- read_csv(file_path, show_col_types = FALSE)

  if (!(lyap_col_name %in% names(df))) {
    stop(paste("Error: Expected column", lyap_col_name, "not found in file:", file_name))
  }
    
  df <- df %>%
    # Filter out LLE failures (where Python used the sentinel value)
    filter(.data[[lyap_col_name]] != LYE_FAILURE_SENTINEL) %>%
    mutate(
      Subject = as.factor(Subject),
      Block = as.factor(Block),
      Trial = as.factor(Trial),
      Difficulty = as.factor(Sequence),
      Trial_RT = trial.RTS, 
      COM_acc = as.character(COM_acc)
    )
  
    df <- df %>%
      arrange(Subject, Block, Trial) %>%
      group_by(Subject) %>%
      mutate(Total_Trial = row_number()) %>% 
      ungroup() %>%
      mutate(log_RT = log(Trial_RT)) %>%
      mutate(
        Phase = case_when(
          Block %in% 1:3 ~ "Learning",
          Block == 4 ~ "Retention",
          Block == 5 ~ "Transfer",
          TRUE ~ "Other" # Safety net
          ),
        Phase = factor(Phase, levels = c("Learning", "Retention", "Transfer"))) %>%
      select(Subject, Block, Phase, Trial, Difficulty, Trial_RT, log_RT, Trial_Acc, COM_acc, all_of(lyap_col_name), Total_Trial)
}

# List to hold data for all participants
all_prep_data_list <- list()
all_exec_data_list <- list()

for (pid in PARTICIPANT_IDS) {
  # Temporarily set the current PID for the function 
  current_prep <- load_lyap_data(pid, "prep")
  current_exec <- load_lyap_data(pid, "exec")
  
  all_prep_data_list[[as.character(pid)]] <- current_prep
  all_exec_data_list[[as.character(pid)]] <- current_exec
}

# Combine all data frames
data_prep_combined <- bind_rows(all_prep_data_list)
data_exec_combined <- bind_rows(all_exec_data_list)

# 2. Prepare Execution Data 
data_exec_renamed <- data_exec_combined %>%
    select(Subject, Block, Trial, Lyap_exec, COM_acc) %>%
    rename(COM_acc_Exec = COM_acc)

# 3. Prepare Preparation Data 
data_prep_renamed <- data_prep_combined %>%
    select(Subject, Block, Trial, Phase, Difficulty, Total_Trial, 
           Trial_RT, log_RT, Trial_Acc, Lyap_prep, COM_acc) %>%
    rename(COM_acc_Prep = COM_acc)

# 4. Final Join and Standardization
data_all_combined <- data_prep_renamed %>%
    left_join(data_exec_renamed, by = c("Subject", "Block", "Trial")) %>%
    mutate(PE = (100 - Trial_Acc)/100)

sd_rt_global <- sd(data_all_combined$Trial_RT, na.rm = TRUE)
sd_pe_global <- sd(data_all_combined$PE, na.rm = TRUE)
LISAS_weight <- sd_rt_global / sd_pe_global

data_all_combined <- data_all_combined %>%
  mutate(LISAS = Trial_RT + (PE *  LISAS_weight)) 

data_all_learning <- data_all_combined %>% filter(Phase == "Learning")
data_all_testing <- data_all_combined %>% filter(Phase %in% c("Retention", "Transfer"))
data_all_retention <- data_all_combined %>% filter(Phase == "Retention")
data_all_transfer <- data_all_combined %>% filter(Phase == "Transfer")
data_all_lr <- data_all_combined %>% filter(Phase %in% c("Learning", "Retention"))

model_LISAS_learn <- lmer(LISAS ~ Difficulty + (1 | Subject), data = data_all_learning)
summary(model_LISAS_learn)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: LISAS ~ Difficulty + (1 | Subject)
##    Data: data_all_learning
## 
## REML criterion at convergence: 3173.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.0055 -0.6145 -0.3021  0.3345 10.6325 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Subject  (Intercept) 0.0843   0.2904  
##  Residual             0.2067   0.4546  
## Number of obs: 2449, groups:  Subject, 18
## 
## Fixed effects:
##               Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)  6.642e-01  7.020e-02 1.816e+01   9.462 1.91e-08 ***
## Difficulty12 8.724e-02  2.213e-02 2.429e+03   3.942 8.30e-05 ***
## Difficulty18 1.610e-01  2.272e-02 2.429e+03   7.084 1.83e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Dffc12
## Difficlty12 -0.157       
## Difficlty18 -0.153  0.486

anova(model_LISAS_learn)

## Type III Analysis of Variance Table with Satterthwaite's method
##            Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
## Difficulty 10.438  5.2188     2  2429  25.253 1.398e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

emmeans(model_LISAS_learn, specs = ~ Difficulty, cov.reduce = range)

##  Difficulty emmean     SE   df lower.CL upper.CL
##  6           0.664 0.0702 18.2    0.517    0.812
##  12          0.751 0.0702 18.2    0.604    0.899
##  18          0.825 0.0704 18.4    0.678    0.973
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95

model_LISAS_test <- lmer(LISAS ~ Difficulty * Block + (1 | Subject), data = data_all_testing)
summary(model_LISAS_test)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: LISAS ~ Difficulty * Block + (1 | Subject)
##    Data: data_all_testing
## 
## REML criterion at convergence: 1793.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.9505 -0.6284 -0.2779  0.3787  5.7532 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Subject  (Intercept) 0.05899  0.2429  
##  Residual             0.16219  0.4027  
## Number of obs: 1674, groups:  Subject, 18
## 
## Fixed effects:
##                      Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)         5.578e-01  6.224e-02 2.241e+01   8.963 7.27e-09 ***
## Difficulty12        5.816e-02  3.435e-02 1.651e+03   1.693  0.09065 .  
## Difficulty18        1.508e-01  3.426e-02 1.651e+03   4.400 1.15e-05 ***
## Block5              1.869e-01  3.417e-02 1.651e+03   5.468 5.24e-08 ***
## Difficulty12:Block5 1.408e-01  4.828e-02 1.651e+03   2.916  0.00359 ** 
## Difficulty18:Block5 1.129e-01  4.824e-02 1.651e+03   2.339  0.01944 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Dffc12 Dffc18 Block5 D12:B5
## Difficlty12 -0.279                            
## Difficlty18 -0.280  0.507                     
## Block5      -0.281  0.508  0.510              
## Dffclt12:B5  0.199 -0.712 -0.361 -0.708       
## Dffclt18:B5  0.199 -0.360 -0.710 -0.709  0.502

anova(model_LISAS_test)

## Type III Analysis of Variance Table with Satterthwaite's method
##                  Sum Sq Mean Sq NumDF DenDF  F value    Pr(>F)    
## Difficulty       12.196  6.0978     2  1651  37.5960 < 2.2e-16 ***
## Block            30.820 30.8198     1  1651 190.0199 < 2.2e-16 ***
## Difficulty:Block  1.546  0.7730     2  1651   4.7661  0.008631 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

emmeans(model_LISAS_test, specs = ~ Difficulty, cov.reduce = range)

## NOTE: Results may be misleading due to involvement in interactions

##  Difficulty emmean     SE df lower.CL upper.CL
##  6           0.651 0.0597 19    0.526    0.776
##  12          0.780 0.0597 19    0.655    0.905
##  18          0.858 0.0597 19    0.733    0.983
## 
## Results are averaged over the levels of: Block 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95

plot_data_LISAS_raw <- data_all_combined %>%
  group_by(Block, Difficulty) %>%
  summarise(
    Num_Trial = as.numeric(Trial),
    LISAS = as.numeric(LISAS),
    mean_LISAS = mean(LISAS, na.rm = TRUE),
    sd_LISAS = sd(LISAS, na.rm = TRUE),
    n = n(),
    se_LISAS = sd_LISAS / sqrt(n)
  )

## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
##   always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## `summarise()` has grouped output by 'Block', 'Difficulty'. You can override
## using the `.groups` argument.

block_names <- c(
  "1" = "Block 1: Learning",
  "2" = "Block 2: Learning",
  "3" = "Block 3: Learning",
  "4" = "Block 4: Retention",
  "5" = "Block 5: Transfer"
)

p_LISAS_raw <- ggplot(plot_data_LISAS_raw, aes(x = Num_Trial, y = LISAS, color = Difficulty, fill = Difficulty)) +
  geom_boxplot(alpha = 0.7, outlier.shape = NA) +
  geom_jitter(width = 0.2, alpha = 0.1) +
  facet_wrap(~ Block, labeller = as_labeller(block_names)) +
  scale_x_continuous(breaks = c(12, 24, 36, 48)) +
  labs(
    title = "Performance Across Learning, Retention, and Transfer Blocks",
    x = "Trial Number",
    y = "LISAS"
    #color = "Difficulty"
  ) +
  theme_classic() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
    axis.title = element_text(size = 14),
    legend.position = "bottom"
  )




plot_data_LISAS_block <- data_all_combined %>%
  group_by(Block) %>%
  summarise(
    mean_LISAS = mean(LISAS, na.rm = TRUE),
    sd_LISAS = sd(LISAS, na.rm = TRUE),
    n = n(),
    se_LISAS = sd_LISAS / sqrt(n)
  )

p_LISAS_block <- ggplot(plot_data_LISAS_block, aes(x = as.factor(Block), y =  mean_LISAS)) +
  geom_vline(xintercept = 3.5, color = "red", linewidth = 1) +
  geom_errorbar(aes(ymin =  mean_LISAS - se_LISAS, 
                    ymax =  mean_LISAS + se_LISAS), 
                width = 0.2, color = "darkblue", linewidth = 0.8) +
  geom_point(size = 4, color = "darkblue") +
  labs(
    title = "LISAS ~ Block",
    x = "Block",
    y = "Mean LISAS"
  ) +
  theme_classic() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 16, face = "plain"),
    axis.title = element_text(size = 14),
    axis.text = element_text(size = 12)
  )

p_LISAS_raw

p_LISAS_block

Individual learning curves

# 1. Prepare Data
# Filter data to only include the Learning Phase (Blocks 1-3)
# Ensure Total_Trial is numeric and Subject and Difficulty are factors
data_learning_individual <- data_all_learning %>%
  mutate(
    Trial = as.numeric(Trial),
    Total_Trial = as.numeric(Total_Trial),
    Subject = as.factor(Subject),
    Difficulty = as.factor(Difficulty)
  )

# 2. Generate Individual Plots using a Loop or facet_wrap
p_individual_exec_learning_curves <- ggplot(data_learning_individual, 
       aes(x = Total_Trial, y = Lyap_exec, color = Difficulty)) +
  geom_point(alpha = 0.4, size = 1) +
  geom_smooth(method = "gam", 
              formula = y ~ s(x), 
              se = FALSE, 
              linewidth = 0.8) +
  facet_wrap(~ Subject, scales = "free_y") +
  scale_color_manual(values = c("6" = "red", "12" = "green", "18" = "blue"),
                     name = "Sequence\nDifficulty") +
  labs(
    title = "Individual Learning Curves: Execution Stability",
    subtitle = "Progression across Learning Trials (Blocks 1-3) by Sequence Difficulty",
    x = "Total Trial Number",
    y = expression("Lyapunov Exponent")
  ) +
  theme_bw(base_size = 10) +
  theme(
    strip.text = element_text(face = "bold"), 
    legend.position = "bottom"
  )
print(p_individual_exec_learning_curves)

p_individual_prep_learning_curves <- ggplot(data_learning_individual, 
       aes(x = Total_Trial, y = Lyap_prep, color = Difficulty)) +
  geom_point(alpha = 0.4, size = 1) +
  geom_smooth(method = "gam", 
              formula = y ~ s(x), 
              se = FALSE,
              linewidth = 0.8) +
  facet_wrap(~ Subject, scales = "free_y") + 
  scale_color_manual(values = c("6" = "red", "12" = "green", "18" = "blue"),
                     name = "Sequence\nDifficulty") +
  labs(
    title = "Individual Learning Curves: Preparation Stability",
    subtitle = "Progression across Learning Trials (Blocks 1-3) by Sequence Difficulty",
    x = "Total Trial Number",
    y = expression("Lyapunov Exponent" )
  ) +
  theme_bw(base_size = 10) +
  theme(
    strip.text = element_text(face = "bold"), 
    legend.position = "bottom"
  )
print(p_individual_prep_learning_curves)

Learning Phase

print("Difficulty effect on Preparation (Learning)")

## [1] "Difficulty effect on Preparation (Learning)"

model_lyap_prep_learning <- lmer(Lyap_prep ~ Difficulty * Block + (1 | Subject), data = data_all_learning)

## fixed-effect model matrix is rank deficient so dropping 6 columns / coefficients

summary(model_lyap_prep_learning)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Lyap_prep ~ Difficulty * Block + (1 | Subject)
##    Data: data_all_learning
## 
## REML criterion at convergence: -13629.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.8554 -0.5196 -0.0662  0.5434  4.8712 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  Subject  (Intercept) 1.049e-05 0.003239
##  Residual             2.177e-04 0.014755
## Number of obs: 2449, groups:  Subject, 18
## 
## Fixed effects:
##               Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)  4.310e-02  9.167e-04 2.658e+01  47.018  < 2e-16 ***
## Difficulty12 2.372e-03  7.183e-04 2.429e+03   3.302 0.000972 ***
## Difficulty18 3.357e-03  7.375e-04 2.429e+03   4.552 5.57e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Dffc12
## Difficlty12 -0.391       
## Difficlty18 -0.381  0.486
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 6 columns / coefficients

anova(model_lyap_prep_learning)

## Missing cells for: Difficulty12:Block1, Difficulty18:Block1, Difficulty6:Block2, Difficulty18:Block2, Difficulty6:Block3, Difficulty12:Block3.  
## Interpret type III hypotheses with care.

## Type III Analysis of Variance Table with Satterthwaite's method
##                     Sum Sq   Mean Sq NumDF  DenDF F value    Pr(>F)    
## Difficulty       0.0048496 0.0024248     2 2429.1  11.137 1.531e-05 ***
## Block                                                                  
## Difficulty:Block                                                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

print("Difficulty effect on Execution (Learning)")

## [1] "Difficulty effect on Execution (Learning)"

model_lyap_exec_learning <- lmer(Lyap_exec ~ Difficulty * Block + (1 | Subject), data = data_all_learning)

## fixed-effect model matrix is rank deficient so dropping 6 columns / coefficients

summary(model_lyap_exec_learning)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Lyap_exec ~ Difficulty * Block + (1 | Subject)
##    Data: data_all_learning
## 
## REML criterion at convergence: -13463.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -6.1935 -0.5928 -0.0478  0.5799  3.3795 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  Subject  (Intercept) 2.908e-05 0.005392
##  Residual             2.316e-04 0.015219
## Number of obs: 2449, groups:  Subject, 18
## 
## Fixed effects:
##                Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)   5.383e-02  1.375e-03  2.075e+01  39.161   <2e-16 ***
## Difficulty12  7.528e-04  7.409e-04  2.429e+03   1.016    0.310    
## Difficulty18 -3.750e-05  7.606e-04  2.429e+03  -0.049    0.961    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Dffc12
## Difficlty12 -0.269       
## Difficlty18 -0.262  0.486
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 6 columns / coefficients

anova(model_lyap_exec_learning)

## Missing cells for: Difficulty12:Block1, Difficulty18:Block1, Difficulty6:Block2, Difficulty18:Block2, Difficulty6:Block3, Difficulty12:Block3.  
## Interpret type III hypotheses with care.

## Type III Analysis of Variance Table with Satterthwaite's method
##                      Sum Sq   Mean Sq NumDF DenDF F value Pr(>F)
## Difficulty       0.00032879 0.0001644     2  2429  0.7098 0.4919
## Block                                                           
## Difficulty:Block

emm_exec_learning <- emmeans(model_lyap_exec_learning, specs = ~ Difficulty, cov.reduce = range)

## NOTE: A nesting structure was detected in the fitted model:
##     Difficulty %in% Block

print(emm_exec_learning)

##  Difficulty Block emmean      SE   df lower.CL upper.CL
##  6          1     0.0538 0.00137 20.8   0.0510   0.0567
##  12         2     0.0546 0.00137 20.8   0.0517   0.0574
##  18         3     0.0538 0.00139 21.4   0.0509   0.0567
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95

emm_prep_learning <- emmeans(model_lyap_prep_learning, specs = ~ Difficulty, cov.reduce = range)

## NOTE: A nesting structure was detected in the fitted model:
##     Difficulty %in% Block

print(emm_prep_learning)

##  Difficulty Block emmean       SE   df lower.CL upper.CL
##  6          1     0.0431 0.000917 26.6   0.0412   0.0450
##  12         2     0.0455 0.000917 26.6   0.0436   0.0474
##  18         3     0.0465 0.000932 28.4   0.0446   0.0484
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95

learn_prep_emms <- emmeans(model_lyap_prep_learning, ~ Difficulty)

## NOTE: A nesting structure was detected in the fitted model:
##     Difficulty %in% Block

learn_prep_df <- as.data.frame(learn_prep_emms)

ggplot(learn_prep_df, aes(x = Difficulty, y = emmean, group = 1)) +
  geom_line(size = 1, color = "darkblue") +
  geom_point(size = 4, color = "darkblue") +
  geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = 0.1, size = 1, color = "black") +
  labs(title = "Preparation Stability Across Sequence Difficulty",
       x = "Sequence Length (Difficulty)",
       y = "Estimated Marginal Mean (sLyE Prep)",
       caption = "Error bars represent 95% Confidence Intervals") +
  theme_classic() +
  theme(text = element_text(size = 14))

## 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.

learn_exec_emms <- emmeans(model_lyap_exec_learning, ~ Difficulty)

## NOTE: A nesting structure was detected in the fitted model:
##     Difficulty %in% Block

learn_exec_df <- as.data.frame(learn_exec_emms)

ggplot(learn_exec_df, aes(x = Difficulty, y = emmean, group = 1)) +
  geom_line(size = 1, color = "darkred") +
  geom_point(size = 4, color = "darkred") +
  geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = 0.1, size = 1, color = "black") +
  labs(title = "Execution Stability Across Sequence Difficulty",
       x = "Sequence Length (Difficulty)",
       y = "Estimated Marginal Mean (sLyE Exec)",
       caption = "Error bars represent 95% Confidence Intervals") +
  theme_classic() +
  theme(text = element_text(size = 14))

Test Phase

model_prep_test <- lmer(Lyap_prep ~ Block * Difficulty + (1 | Subject),
                        data = data_all_testing)
summary(model_prep_test)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Lyap_prep ~ Block * Difficulty + (1 | Subject)
##    Data: data_all_testing
## 
## REML criterion at convergence: -9177.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.4618 -0.5903 -0.0611  0.5932  4.1403 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  Subject  (Intercept) 1.921e-05 0.004382
##  Residual             2.289e-04 0.015128
## Number of obs: 1674, groups:  Subject, 18
## 
## Fixed effects:
##                       Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)          4.608e-02  1.382e-03  4.258e+01  33.350  < 2e-16 ***
## Block5               7.646e-04  1.284e-03  1.651e+03   0.596  0.55152    
## Difficulty12         3.522e-03  1.290e-03  1.651e+03   2.730  0.00641 ** 
## Difficulty18         2.818e-04  1.287e-03  1.651e+03   0.219  0.82670    
## Block5:Difficulty12 -8.599e-04  1.813e-03  1.651e+03  -0.474  0.63542    
## Block5:Difficulty18  7.956e-05  1.812e-03  1.651e+03   0.044  0.96498    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Block5 Dffc12 Dffc18 B5:D12
## Block5      -0.475                            
## Difficlty12 -0.472  0.508                     
## Difficlty18 -0.473  0.510  0.507              
## Blck5:Dff12  0.336 -0.708 -0.712 -0.361       
## Blck5:Dff18  0.336 -0.709 -0.360 -0.710  0.501

anova(model_prep_test)

## Type III Analysis of Variance Table with Satterthwaite's method
##                     Sum Sq    Mean Sq NumDF  DenDF F value    Pr(>F)    
## Block            0.0001065 0.00010646     1 1651.1  0.4652 0.4953172    
## Difficulty       0.0032232 0.00161160     2 1651.1  7.0420 0.0009009 ***
## Block:Difficulty 0.0000757 0.00003787     2 1651.1  0.1655 0.8475233    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

model_exec_test <- lmer(Lyap_exec ~ Block * Difficulty + (1 | Subject),
                        data = data_all_testing)
summary(model_exec_test)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Lyap_exec ~ Block * Difficulty + (1 | Subject)
##    Data: data_all_testing
## 
## REML criterion at convergence: -9116.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.1105 -0.6538 -0.0588  0.6022  3.6113 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  Subject  (Intercept) 2.279e-05 0.004774
##  Residual             2.371e-04 0.015398
## Number of obs: 1674, groups:  Subject, 18
## 
## Fixed effects:
##                       Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)          5.425e-02  1.462e-03  3.912e+01  37.098  < 2e-16 ***
## Block5               3.547e-03  1.307e-03  1.651e+03   2.714  0.00671 ** 
## Difficulty12         3.428e-04  1.313e-03  1.651e+03   0.261  0.79411    
## Difficulty18         1.267e-03  1.310e-03  1.651e+03   0.968  0.33340    
## Block5:Difficulty12  7.864e-04  1.846e-03  1.651e+03   0.426  0.67013    
## Block5:Difficulty18 -1.998e-03  1.844e-03  1.651e+03  -1.083  0.27890    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Block5 Dffc12 Dffc18 B5:D12
## Block5      -0.457                            
## Difficlty12 -0.454  0.508                     
## Difficlty18 -0.455  0.510  0.507              
## Blck5:Dff12  0.323 -0.708 -0.712 -0.361       
## Blck5:Dff18  0.323 -0.709 -0.360 -0.710  0.501

anova(model_exec_test)

## Type III Analysis of Variance Table with Satterthwaite's method
##                     Sum Sq   Mean Sq NumDF  DenDF F value    Pr(>F)    
## Block            0.0041328 0.0041328     1 1651.1 17.4304 3.136e-05 ***
## Difficulty       0.0001545 0.0000773     2 1651.1  0.3259    0.7219    
## Block:Difficulty 0.0005758 0.0002879     2 1651.1  1.2142    0.2972    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

test_prep_emms <- emmeans(model_prep_test, ~ Block * Difficulty)
test_prep_df <- as.data.frame(test_prep_emms)

block_names <- c(
  "1" = "Block 1: Learning",
  "2" = "Block 2: Learning",
  "3" = "Block 3: Learning",
  "4" = "Block 4: Retention",
  "5" = "Block 5: Transfer"
)

ggplot(test_prep_df, aes(x = Difficulty, y = emmean, group = 1, color = Difficulty)) +
  geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = 0.1, size = 1, color = "black") +
  geom_line(size = 1, color = "gray50") + 
  geom_point(size = 4) +
  facet_wrap(~Block, labeller = as_labeller(block_names)) +
  scale_color_brewer(palette = "Set1") +
  labs(title = "Preparation Stability Across Difficulties: Retention vs. Transfer",
       subtitle = "Model-estimated means by Difficulty",
       x = "Sequence Length (Difficulty)",
       y = "Estimated Marginal Mean (sLyE Prep)",
       caption = "Error bars represent 95% Confidence Intervals") +
  theme_classic() +
  theme(text = element_text(size = 14),
        legend.position = "right")

test_exec_emms <- emmeans(model_exec_test, ~ Block * Difficulty)
test_exec_df <- as.data.frame(test_exec_emms)

ggplot(test_exec_df, aes(x = Difficulty, y = emmean, group = 1, color = Difficulty)) +
  geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = 0.1, size = 1, color = "black") +
  geom_line(size = 1, color = "gray50") + 
  geom_point(size = 4) +
  facet_wrap(~Block, labeller = as_labeller(block_names)) +
  scale_color_brewer(palette = "Set1") +
  labs(title = "Execution Stability Across Difficulties: Retention vs. Transfer",
       subtitle = "Model-estimated means by Difficulty",
       x = "Sequence Length (Difficulty)",
       y = "Estimated Marginal Mean (sLyE Exec)",
       caption = "Error bars represent 95% Confidence Intervals") +
  theme_classic() +
  theme(text = element_text(size = 14),
        legend.position = "right")

LISAS

model_lisas_lr <- lmer(LISAS ~ Lyap_prep * Lyap_exec * Difficulty * Block + (1 | Subject), data = data_all_lr)

## fixed-effect model matrix is rank deficient so dropping 24 columns / coefficients

summary(model_lisas_lr)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: LISAS ~ Lyap_prep * Lyap_exec * Difficulty * Block + (1 | Subject)
##    Data: data_all_lr
## 
## REML criterion at convergence: 3812.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.9446 -0.6052 -0.2890  0.2878 11.0322 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Subject  (Intercept) 0.07624  0.2761  
##  Residual             0.18766  0.4332  
## Number of obs: 3280, groups:  Subject, 18
## 
## Fixed effects:
##                                   Estimate Std. Error        df t value
## (Intercept)                         0.8408     0.1715  674.3789   4.903
## Lyap_prep                          -0.9578     3.6019 3240.1200  -0.266
## Lyap_exec                          -5.5334     2.8056 3240.5581  -1.972
## Difficulty12                        0.6002     0.4304 3239.6171   1.394
## Difficulty18                        0.8613     0.4851 3239.2934   1.776
## Block2                             -0.9523     0.4952 3239.5953  -1.923
## Block3                             -1.0210     0.5377 3239.2426  -1.899
## Block4                             -0.8251     0.3517 3239.4295  -2.346
## Lyap_prep:Lyap_exec                69.9553    63.6511 3240.0758   1.099
## Lyap_prep:Difficulty12             -8.6969     8.3860 3239.6431  -1.037
## Lyap_prep:Difficulty18            -13.4878    10.1243 3239.3872  -1.332
## Lyap_exec:Difficulty12            -10.1262     7.5201 3239.6052  -1.347
## Lyap_exec:Difficulty18            -14.2107     8.7546 3239.2425  -1.623
## Lyap_prep:Block2                   16.9582    10.0023 3239.7155   1.695
## Lyap_prep:Block3                   18.0233    11.2822 3239.3762   1.598
## Lyap_prep:Block4                   10.2336     7.1037 3239.6161   1.441
## Lyap_exec:Block2                   20.1954     8.6831 3239.5567   2.326
## Lyap_exec:Block3                   21.9466     9.6750 3239.2148   2.268
## Lyap_exec:Block4                   15.6299     6.3983 3239.4637   2.443
## Lyap_prep:Lyap_exec:Difficulty12  164.2475   148.2606 3239.6408   1.108
## Lyap_prep:Lyap_exec:Difficulty18  273.1538   182.0302 3239.3216   1.501
## Lyap_prep:Lyap_exec:Block2       -360.9952   176.6519 3239.6684  -2.044
## Lyap_prep:Lyap_exec:Block3       -399.5020   202.3394 3239.3386  -1.974
## Lyap_prep:Lyap_exec:Block4       -244.9103   130.8697 3239.6462  -1.871
##                                  Pr(>|t|)    
## (Intercept)                      1.18e-06 ***
## Lyap_prep                          0.7903    
## Lyap_exec                          0.0487 *  
## Difficulty12                       0.1633    
## Difficulty18                       0.0759 .  
## Block2                             0.0546 .  
## Block3                             0.0577 .  
## Block4                             0.0190 *  
## Lyap_prep:Lyap_exec                0.2718    
## Lyap_prep:Difficulty12             0.2998    
## Lyap_prep:Difficulty18             0.1829    
## Lyap_exec:Difficulty12             0.1782    
## Lyap_exec:Difficulty18             0.1046    
## Lyap_prep:Block2                   0.0901 .  
## Lyap_prep:Block3                   0.1103    
## Lyap_prep:Block4                   0.1498    
## Lyap_exec:Block2                   0.0201 *  
## Lyap_exec:Block3                   0.0234 *  
## Lyap_exec:Block4                   0.0146 *  
## Lyap_prep:Lyap_exec:Difficulty12   0.2680    
## Lyap_prep:Lyap_exec:Difficulty18   0.1336    
## Lyap_prep:Lyap_exec:Block2         0.0411 *  
## Lyap_prep:Lyap_exec:Block3         0.0484 *  
## Lyap_prep:Lyap_exec:Block4         0.0614 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 
## Correlation matrix not shown by default, as p = 24 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it

## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 24 columns / coefficients

anova(model_lisas_lr)

## Missing cells for: Difficulty12:Block1, Difficulty18:Block1, Difficulty6:Block2, Difficulty18:Block2, Difficulty6:Block3, Difficulty12:Block3, Lyap_prep:Difficulty12:Block1, Lyap_prep:Difficulty18:Block1, Lyap_prep:Difficulty6:Block2, Lyap_prep:Difficulty18:Block2, Lyap_prep:Difficulty6:Block3, Lyap_prep:Difficulty12:Block3, Lyap_exec:Difficulty12:Block1, Lyap_exec:Difficulty18:Block1, Lyap_exec:Difficulty6:Block2, Lyap_exec:Difficulty18:Block2, Lyap_exec:Difficulty6:Block3, Lyap_exec:Difficulty12:Block3, Lyap_prep:Lyap_exec:Difficulty12:Block1, Lyap_prep:Lyap_exec:Difficulty18:Block1, Lyap_prep:Lyap_exec:Difficulty6:Block2, Lyap_prep:Lyap_exec:Difficulty18:Block2, Lyap_prep:Lyap_exec:Difficulty6:Block3, Lyap_prep:Lyap_exec:Difficulty12:Block3.  
## Interpret type III hypotheses with care.

## Type III Analysis of Variance Table with Satterthwaite's method
##                                       Sum Sq Mean Sq NumDF  DenDF F value
## Lyap_prep                            0.39851 0.39851     1 3240.8  2.1236
## Lyap_exec                            0.05929 0.05929     1 3243.0  0.3160
## Difficulty                           0.66581 0.33290     2 3239.5  1.7740
## Block                                1.10111 0.36704     3 3239.4  1.9559
## Lyap_prep:Lyap_exec                  0.19287 0.19287     1 3241.2  1.0278
## Lyap_prep:Difficulty                 0.37996 0.18998     2 3239.5  1.0124
## Lyap_exec:Difficulty                 0.56492 0.28246     2 3239.4  1.5052
## Lyap_prep:Block                      0.70892 0.23631     3 3239.4  1.2593
## Lyap_exec:Block                      1.37609 0.45870     3 3239.4  2.4443
## Lyap_prep:Lyap_exec:Difficulty       0.45815 0.22907     2 3239.5  1.2207
## Lyap_prep:Lyap_exec:Block            1.02564 0.34188     3 3239.5  1.8218
## Difficulty:Block                                                         
## Lyap_prep:Difficulty:Block                                               
## Lyap_exec:Difficulty:Block                                               
## Lyap_prep:Lyap_exec:Difficulty:Block                                     
##                                       Pr(>F)  
## Lyap_prep                            0.14514  
## Lyap_exec                            0.57408  
## Difficulty                           0.16982  
## Block                                0.11846  
## Lyap_prep:Lyap_exec                  0.31075  
## Lyap_prep:Difficulty                 0.36346  
## Lyap_exec:Difficulty                 0.22213  
## Lyap_prep:Block                      0.28667  
## Lyap_exec:Block                      0.06221 .
## Lyap_prep:Lyap_exec:Difficulty       0.29516  
## Lyap_prep:Lyap_exec:Block            0.14094  
## Difficulty:Block                              
## Lyap_prep:Difficulty:Block                    
## Lyap_exec:Difficulty:Block                    
## Lyap_prep:Lyap_exec:Difficulty:Block          
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

model_lisas_transfer <- lmer(LISAS ~ Lyap_prep * Lyap_exec * Difficulty + (1 | Subject), data = data_all_transfer)
summary(model_lisas_transfer)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: LISAS ~ Lyap_prep * Lyap_exec * Difficulty + (1 | Subject)
##    Data: data_all_transfer
## 
## REML criterion at convergence: 960.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.9502 -0.6701 -0.2827  0.5317  4.9577 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Subject  (Intercept) 0.07638  0.2764  
##  Residual             0.18152  0.4261  
## Number of obs: 843, groups:  Subject, 18
## 
## Fixed effects:
##                                   Estimate Std. Error        df t value
## (Intercept)                        0.61018    0.32722 809.39327   1.865
## Lyap_prep                          0.79771    6.63813 816.37953   0.120
## Lyap_exec                          3.21092    5.33416 816.58576   0.602
## Difficulty12                       0.58154    0.42576 815.67351   1.366
## Difficulty18                       0.24199    0.41588 817.14757   0.582
## Lyap_prep:Lyap_exec              -33.01318  112.48191 815.97206  -0.293
## Lyap_prep:Difficulty12           -10.74872    8.59110 815.60562  -1.251
## Lyap_prep:Difficulty18             0.08581    8.79569 816.31971   0.010
## Lyap_exec:Difficulty12            -8.90200    7.12980 815.83923  -1.249
## Lyap_exec:Difficulty18            -0.80852    6.83190 816.99994  -0.118
## Lyap_prep:Lyap_exec:Difficulty12 231.74975  145.74994 815.69611   1.590
## Lyap_prep:Lyap_exec:Difficulty18  25.00347  146.08705 816.19676   0.171
##                                  Pr(>|t|)  
## (Intercept)                        0.0626 .
## Lyap_prep                          0.9044  
## Lyap_exec                          0.5474  
## Difficulty12                       0.1723  
## Difficulty18                       0.5608  
## Lyap_prep:Lyap_exec                0.7692  
## Lyap_prep:Difficulty12             0.2112  
## Lyap_prep:Difficulty18             0.9922  
## Lyap_exec:Difficulty12             0.2122  
## Lyap_exec:Difficulty18             0.9058  
## Lyap_prep:Lyap_exec:Difficulty12   0.1122  
## Lyap_prep:Lyap_exec:Difficulty18   0.8641  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Lyp_pr Lyp_xc Dffc12 Dffc18 Ly_:L_ Lyp_p:D12 Lyp_p:D18
## Lyap_prep   -0.937                                                       
## Lyap_exec   -0.944  0.930                                                
## Difficlty12 -0.731  0.715  0.720                                         
## Difficlty18 -0.753  0.734  0.742  0.575                                  
## Lyp_prp:Ly_  0.894 -0.964 -0.956 -0.684 -0.703                           
## Lyp_prp:D12  0.720 -0.769 -0.716 -0.951 -0.564  0.743                    
## Lyp_prp:D18  0.706 -0.752 -0.703 -0.540 -0.951  0.727  0.578             
## Lyp_xc:Df12  0.702 -0.693 -0.745 -0.961 -0.555  0.713  0.921     0.526   
## Lyp_xc:Df18  0.735 -0.725 -0.780 -0.563 -0.958  0.746  0.558     0.916   
## Lyp_:L_:D12 -0.690  0.743  0.737  0.909  0.543 -0.771 -0.962    -0.561   
## Lyp_:L_:D18 -0.689  0.742  0.737  0.528  0.907 -0.770 -0.571    -0.961   
##             Lyp_x:D12 Lyp_x:D18 L_:L_:D12
## Lyap_prep                                
## Lyap_exec                                
## Difficlty12                              
## Difficlty18                              
## Lyp_prp:Ly_                              
## Lyp_prp:D12                              
## Lyp_prp:D18                              
## Lyp_xc:Df12                              
## Lyp_xc:Df18  0.584                       
## Lyp_:L_:D12 -0.953    -0.576             
## Lyp_:L_:D18 -0.552    -0.949     0.594

anova(model_lisas_transfer)

## Type III Analysis of Variance Table with Satterthwaite's method
##                                 Sum Sq  Mean Sq NumDF  DenDF F value Pr(>F)
## Lyap_prep                      0.11453 0.114534     1 816.83  0.6310 0.4272
## Lyap_exec                      0.00002 0.000016     1 817.51  0.0001 0.9926
## Difficulty                     0.34993 0.174963     2 816.49  0.9639 0.3819
## Lyap_prep:Lyap_exec            0.15061 0.150610     1 816.08  0.8297 0.3626
## Lyap_prep:Difficulty           0.43088 0.215442     2 816.10  1.1869 0.3057
## Lyap_exec:Difficulty           0.38583 0.192917     2 816.35  1.0628 0.3460
## Lyap_prep:Lyap_exec:Difficulty 0.62667 0.313337     2 815.95  1.7262 0.1786

emmeans(model_lisas_transfer, ~ Difficulty)

## NOTE: Results may be misleading due to involvement in interactions

##  Difficulty emmean   SE   df lower.CL upper.CL
##  6           0.743 0.07 20.4    0.597    0.889
##  12          0.937 0.07 20.5    0.791    1.083
##  18          1.011 0.07 20.5    0.866    1.157
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95

lr_lisas_emms <- emmeans(model_lisas_lr, ~ Block * Difficulty)

## Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
## To enable adjustments, add the argument 'pbkrtest.limit = 3280' (or larger)
## [or, globally, 'set emm_options(pbkrtest.limit = 3280)' 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 = 3280' (or larger)
## [or, globally, 'set emm_options(lmerTest.limit = 3280)' or larger];
## but be warned that this may result in large computation time and memory use.

## NOTE: Results may be misleading due to involvement in interactions

lr_lisas_df <- as.data.frame(lr_lisas_emms)

lr_lisas_df$Block <- factor(lr_lisas_df$Block)
lr_lisas_df$Difficulty <- factor(lr_lisas_df$Difficulty)

pd <- position_dodge(0.5)

# 2. Create the combined plot
ggplot(lr_lisas_df, aes(x = Block, y = emmean, group = Difficulty, color = Difficulty)) +
  geom_errorbar(aes(ymin = emmean - SE, ymax = emmean + SE), width = 0.2, position = pd, alpha = 0.7) +
  geom_point(position = pd, size = 3) +
  geom_vline(xintercept = 3.5, linetype = "dashed", color = "gray40") +
  annotate("text", x = 2, y = max(lr_lisas_df$emmean), label = "Learning Phase", fontface = "italic") +
  annotate("text", x = 4.5, y = max(lr_lisas_df$emmean), label = "Testing Phase", fontface = "italic") +
  scale_color_brewer(palette = "Set1") +
  scale_x_discrete(labels = c("1" = "Learning Block 1", "2" = "Learning Block 2", "3" = "Learning Block 3", "4" = "Retention Block")) +
  labs(title = "Motor Performance (LISAS) from Learning to Retention",
       #subtitle = "Comparison of performance efficiency across the Learning Blocks to the Retention Block",
       x = "Block",
       y = "Mean LISAS (lower is better)") +
       # caption = "Error bars represent ±1 Standard Error. Learning Blocks 1-3: Constant Difficulty; Retention Block: Mixed Difficulty.") +
  theme_classic() +
  theme(text = element_text(size = 14), legend.position = "top")

## Warning: Removed 6 rows containing missing values or values outside the scale range
## (`geom_point()`).

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_text()`).
## Removed 1 row containing missing values or values outside the scale range
## (`geom_text()`).

transfer_lisas_emms <- emmeans(model_lisas_transfer, ~ Difficulty)

## NOTE: Results may be misleading due to involvement in interactions

transfer_lisas_df <- as.data.frame(transfer_lisas_emms)
transfer_lisas_df$Difficulty <- factor(transfer_lisas_df$Difficulty)
  
ggplot(transfer_lisas_df, aes(x = Difficulty, y = emmean, color = Difficulty), group = 1) +
  geom_line(size = 1, color = "darkblue") +
  geom_point(size = 4, color = "darkblue") +
  geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = 0.1, size = 1, color = "black") +
  #scale_x_discrete(labels = c("1" = "Learning Block 1", "2" = "Learning Block 2", "3" = "Learning Block 3", "4" = "Retention")) +
  labs(title = "Motor Performance (LISAS) During the Transfer Block",
       #subtitle = "Comparison of performance efficiency across the Transfer Block",
       x = "Sequence Length (Difficulty)",
       y = "Mean LISAS (lower is better)") +
       # caption = "Error bars represent ±1 Standard Error. Learning Blocks 1-3: Constant Difficulty; Retention Block: Mixed Difficulty.") +
  theme_classic() +
  theme(text = element_text(size = 14), legend.position = "top")

## Warning in fortify(data, ...): Arguments in `...` must be used.
## ✖ Problematic argument:
## • group = 1
## ℹ Did you misspell an argument name?

## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?

sd_p <- sd(data_all_lr$Lyap_prep, na.rm = TRUE)
sd_e <- sd(data_all_lr$Lyap_exec, na.rm = TRUE)

emm_prep <- emmeans(model_lisas_lr, ~ Block + Difficulty, at = list(Lyap_prep = mean(data_all_lr$Lyap_prep) + sd_p))

## Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
## To enable adjustments, add the argument 'pbkrtest.limit = 3280' (or larger)
## [or, globally, 'set emm_options(pbkrtest.limit = 3280)' 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 = 3280' (or larger)
## [or, globally, 'set emm_options(lmerTest.limit = 3280)' or larger];
## but be warned that this may result in large computation time and memory use.

## NOTE: Results may be misleading due to involvement in interactions

emm_exec <- emmeans(model_lisas_lr, ~ Block + Difficulty, at = list(Lyap_exec = mean(data_all_lr$Lyap_exec) + sd_e))

## Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
## To enable adjustments, add the argument 'pbkrtest.limit = 3280' (or larger)
## [or, globally, 'set emm_options(pbkrtest.limit = 3280)' 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 = 3280' (or larger)
## [or, globally, 'set emm_options(lmerTest.limit = 3280)' or larger];
## but be warned that this may result in large computation time and memory use.

## NOTE: Results may be misleading due to involvement in interactions

df_prep <- as.data.frame(emm_prep) %>% mutate(Stability_Type = "Prep sLyE")
df_exec <- as.data.frame(emm_exec) %>% mutate(Stability_Type = "Exec sLyE")
plot_df <- bind_rows(df_prep, df_exec)

ggplot(plot_df, aes(x = Block, y = emmean, group = Stability_Type, color = Stability_Type)) +
  geom_errorbar(aes(ymin = emmean - SE, ymax = emmean + SE), width = 0.2, position = position_dodge(0.4)) +
  geom_point(position = position_dodge(0.4), size = 3) +
  geom_vline(xintercept = 3.5, linetype = "dashed", color = "gray40") +
  facet_wrap(~Difficulty, labeller = labeller(Difficulty = c("6" = "6-Step Difficulty", "12" = "12-Step Difficulty", "18" = "18-Step Difficulty"))) +
  scale_color_manual(values = c("darkred", "darkgreen")) + 
  scale_x_discrete(labels = c("1" = "B1", "2" = "B2", "3" = "B3", "4" = "Ret")) +
  labs(title = "Interaction of sLyE and Difficulty on LISAS Across Learning and Retention Blocks",
       #subtitle = "Comparing how planning vs. execution instability impacts efficiency",
       x = "Block",
       y = "Predicted LISAS",
       color = "Lyapunov Exponent") +
  theme_classic() +
  theme(legend.position = "top", text = element_text(size = 12))

## Warning: Removed 12 rows containing missing values or values outside the scale range
## (`geom_point()`).

sd_p_transfer <- sd(data_all_transfer$Lyap_prep, na.rm = TRUE)
sd_e_transfer <- sd(data_all_transfer$Lyap_exec, na.rm = TRUE)

emm_prep_transfer <- emmeans(model_lisas_transfer, ~ Difficulty, at = list(Lyap_prep = mean(data_all_transfer$Lyap_prep) + sd_p_transfer))

## NOTE: Results may be misleading due to involvement in interactions

emm_exec_transfer <- emmeans(model_lisas_transfer, ~ Difficulty, at = list(Lyap_exec = mean(data_all_transfer$Lyap_exec) + sd_e_transfer))

## NOTE: Results may be misleading due to involvement in interactions

df_prep_transfer <- as.data.frame(emm_prep_transfer) %>% mutate(Stability_Type = "Prep sLyE")
df_exec_transfer <- as.data.frame(emm_exec_transfer) %>% mutate(Stability_Type = "Exec sLyE")
plot_df_transfer <- bind_rows(df_prep_transfer, df_exec_transfer)

ggplot(plot_df_transfer, aes(x = Difficulty, y = emmean, group = Stability_Type, color = Stability_Type)) +
  geom_errorbar(aes(ymin = emmean - SE, ymax = emmean + SE), width = 0.2, position = position_dodge(0.4)) +
  geom_point(position = position_dodge(0.4), size = 3) +
  scale_color_manual(values = c("darkred", "darkgreen")) + 
  labs(title = "Interaction of sLyE and Difficulty on LISAS Across the Transfer Block",
       x = "Sequence Length (Difficulty)",
       y = "Predicted LISAS",
       color = "Lyapunov Exponent") +
  theme_classic() +
  theme(legend.position = "top", text = element_text(size = 12))

Steady states

# 1. Prepare the scaled data
data_for_k <- data_all_combined %>%
  select(Lyap_exec, Lyap_prep, LISAS) %>%
  na.omit() %>%
  scale()

# 2. The Elbow Method (Total Within-Cluster Sum of Squares)
# This looks for the "elbow" where adding more clusters doesn't significantly reduce variance
plot_elbow <- fviz_nbclust(data_for_k, kmeans, method = "wss") +
  labs(title = "Elbow Method: Optimal k for Performance-Stability States",
       subtitle = "Look for the point where the drop flattens out")

# 3. The Silhouette Method (Cluster Quality)
# This measures how similar an object is to its own cluster compared to other clusters
plot_silhouette <- fviz_nbclust(data_for_k, kmeans, method = "silhouette") +
  labs(title = "Silhouette Method: Cluster Consistency",
       subtitle = "Higher values indicate better-defined clusters")

print(plot_elbow) 

print(plot_silhouette)

CLusters

# 1. Define Parameters 
k_optimal <- 4      # Number of steady states (attractors)
set.seed(42)          # For reproducibility

# 2. Prepare Data & Perform K-means Clustering
data_clustering <- data_all_lr %>%
  select(Lyap_exec, Lyap_prep, LISAS) %>%
  na.omit()

# Scale the data (ensures zero variance issues)
data_scaled <- scale(data_clustering)
km_result <- kmeans(data_scaled, centers = k_optimal, nstart = 25) 

# Add cluster assignment back to the data frame
data_clustering$State_ID <- as.factor(km_result$cluster)

# 3. Characterize the Attractors
# This allows us to see which State corresponds to "Best Performance"
attractor_summary <- data_clustering %>%
  group_by(State_ID) %>%
  summarise(
    N = n(),
    mean_Lyap_exec = mean(Lyap_exec),
    mean_Lyap_prep = mean(Lyap_prep),
    mean_LISAS = mean(LISAS)
  ) %>%
  ungroup() %>%
  arrange(State_ID)

print("Learning: Attractor Characterization (k=4) - Sorted by Performance (LISAS)")

## [1] "Learning: Attractor Characterization (k=4) - Sorted by Performance (LISAS)"

print(attractor_summary)

## # A tibble: 4 × 5
##   State_ID     N mean_Lyap_exec mean_Lyap_prep mean_LISAS
##   <fct>    <int>          <dbl>          <dbl>      <dbl>
## 1 1         1162         0.0454         0.0356      0.486
## 2 2          752         0.0731         0.0425      0.596
## 3 3          552         0.0501         0.0457      1.63 
## 4 4          814         0.0523         0.0625      0.541

# 4. Visualization of the Performance-Stability Phase Space
# Mapping LISAS to the size or shape helps visualize how efficiency lives in the stability space
plot_attractor <- ggplot(data_clustering, aes(x = Lyap_exec, y = Lyap_prep, color = State_ID)) + 
  geom_point(alpha = 0.3) +
  # Centroids marked as X
  geom_point(data = as.data.frame(km_result$centers) %>% 
               mutate(State_ID = as.factor(1:k_optimal)), 
             aes(x = Lyap_exec, y = Lyap_prep),
             color = "black", size = 6, shape = 4, stroke = 2) +
  labs(
    title = "Performance-Stability Attractors",
    subtitle = "Clusters defined by Execution Stability, Prep Stability, and LISAS Efficiency",
    x = "Execution Lyapunov Exponent",
    y = "Preparation Lyapunov Exponent"
  ) +
  scale_color_brewer(palette = "Dark2") +
  theme_minimal() +
  theme(plot.title = element_text(face = "bold"))

print(plot_attractor)

# Create the 3D scatter plot
fig_lr <- plot_ly(data_clustering, 
               x = ~Lyap_exec, 
               y = ~Lyap_prep, 
               z = ~LISAS, 
               color = ~State_ID, 
               colors = "Dark2",
               type = 'scatter3d', 
               mode = 'markers',
               marker = list(size = 3, opacity = 0.5)) %>%
  layout(title = "State-Space: Learning & Retention Phase",
         scene = list(xaxis = list(title = 'Exec Stability (sLyE)'),
                      yaxis = list(title = 'Prep Stability (sLyE)'),
                      zaxis = list(title = 'LISAS (Efficiency)')))

fig_lr