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
# Configuration (must match Python script settings)
PARTICIPANT_ID <- 2 # Placeholder for the participant ID 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

Data loading

load_lyap_data <- function(pid, period) {
  # Maps period to folder and file name part
  period_to_folder_map <- list(
    "prep" = "Movement Preparation Period",
    "exec" = "Movement Execution Period"
  )
  period_to_name_map <- list(
    "prep" = "Prep",
    "exec" = "Exec"
  )
  
  folder_name <- period_to_folder_map[[period]]
  name_part <- period_to_name_map[[period]]
  
  # Construct the full path based on the Python saving logic
  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) %>%
    # Filter out LLE failures (where Python used the sentinel value)
    filter(Lyap_Acc != LYE_FAILURE_SENTINEL) %>%
    # Convert necessary columns
    mutate(
      Subject = as.factor(Subject),
      Block = as.factor(Block),
      Trial = as.factor(Trial),
      Sequence = as.factor(Sequence)
    )
    df <- df %>%
      # First, ensure the ordering based on the trial sequence identifiers
      arrange(Subject, Block, Trial) %>%
      # Group by Subject to calculate the running trial number within each participant
      group_by(Subject) %>%
      mutate(Total_Trial = row_number()) %>% # Calculates 1, 2, 3, ... for each subject
      ungroup() %>%
      # Calculate the required log term for the LMEMs
      mutate(log_Total_Trial = log(Total_Trial)) %>%
      mutate(
        Phase = case_when(
          Block %in% 1:3 ~ "Learning",
          Block %in% 4:5 ~ "Testing",
          TRUE ~ "Other" # Safety net
          ),
        Phase = factor(Phase, levels = c("Learning", "Testing"))) %>%
      select(Subject, Block, Phase, Trial, Sequence, trial.RTS, RMS_acc, Lyap_Acc, Total_Trial, log_Total_Trial)
  
  message(paste("Data loaded for", period, ":", nrow(df), "valid rows. Total_Trial and log_Total_Trial successfully calculated."))
  return(df)
}

# Load the datasets
prep_data <- load_lyap_data(PARTICIPANT_ID, "prep")
## Loading data from: /Users/vivi/UT/Thesis/Data/Xsens/Participants/Movement Preparation Period/ID2_Prep_LyE.csv
## Data loaded for prep : 198 valid rows. Total_Trial and log_Total_Trial successfully calculated.
exec_data <- load_lyap_data(PARTICIPANT_ID, "exec")
## Loading data from: /Users/vivi/UT/Thesis/Data/Xsens/Participants/Movement Execution Period/ID2_Exec_LyE.csv
## Data loaded for exec : 198 valid rows. Total_Trial and log_Total_Trial successfully calculated.
# Filter dataframes for specific phases
exec_data_learning <- exec_data %>% filter(Phase == "Learning")
exec_data_testing <- exec_data %>% filter(Phase == "Testing")
prep_data_learning <- prep_data %>% filter(Phase == "Learning")
prep_data_testing <- prep_data %>% filter(Phase == "Testing")

Execution Stability analysis

# All trials (Blocks 1-5)
message("LM for Execution Stability (Lyap_Acc ~ Block * Trial) - Single Subject (B1-5)")
## LM for Execution Stability (Lyap_Acc ~ Block * Trial) - Single Subject (B1-5)
  model_exec_block_trial <- lm(
    Lyap_Acc ~ Block * log_Total_Trial, 
    data = exec_data
  )
  
  message("--- LM Summary for Lyap_Acc ~ Block * Trial (B1-5) ---")
## --- LM Summary for Lyap_Acc ~ Block * Trial (B1-5) ---
  summary(model_exec_block_trial)
## 
## Call:
## lm(formula = Lyap_Acc ~ Block * log_Total_Trial, data = exec_data)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0095302 -0.0016791  0.0001793  0.0017986  0.0064531 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             0.0090365  0.0015014   6.019 9.03e-09 ***
## Block2                  0.0205725  0.0083061   2.477  0.01414 *  
## Block3                  0.0668457  0.0210501   3.176  0.00175 ** 
## Block4                  0.0569317  0.0226317   2.516  0.01272 *  
## Block5                  0.1000726  0.0463573   2.159  0.03214 *  
## log_Total_Trial         0.0009696  0.0005238   1.851  0.06573 .  
## Block2:log_Total_Trial -0.0049048  0.0020591  -2.382  0.01822 *  
## Block3:log_Total_Trial -0.0141975  0.0045788  -3.101  0.00223 ** 
## Block4:log_Total_Trial -0.0118121  0.0045897  -2.574  0.01083 *  
## Block5:log_Total_Trial -0.0194398  0.0089255  -2.178  0.03065 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.002814 on 188 degrees of freedom
## Multiple R-squared:  0.2226, Adjusted R-squared:  0.1854 
## F-statistic: 5.983 on 9 and 188 DF,  p-value: 2.335e-07
# Learning vs Testing blocks
message("3b. Linear Model for Execution Stability - LEARNING (Blocks 1-3) VS TESTING (Blocks 4-5)")
## 3b. Linear Model for Execution Stability - LEARNING (Blocks 1-3) VS TESTING (Blocks 4-5)
# Model compares the mean Lyap_Acc across the block-specific phases, adjusting for sequence.
model_exec_phase_seq <- lm(
  Lyap_Acc ~ Phase * Sequence, 
  data = exec_data
)

message("--- LM Summary for Learning vs Testing Comparison ---")
## --- LM Summary for Learning vs Testing Comparison ---
summary(model_exec_phase_seq)
## 
## Call:
## lm(formula = Lyap_Acc ~ Phase * Sequence, data = exec_data)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0106524 -0.0014421  0.0003411  0.0018104  0.0069897 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              1.169e-02  4.473e-04  26.132  < 2e-16 ***
## PhaseTesting            -1.178e-03  6.849e-04  -1.720   0.0871 .  
## Sequence12               1.800e-03  6.143e-04   2.930   0.0038 ** 
## Sequence18               3.151e-03  6.456e-04   4.881 2.21e-06 ***
## PhaseTesting:Sequence12  1.753e-03  9.671e-04   1.813   0.0715 .  
## PhaseTesting:Sequence18 -7.559e-05  1.012e-03  -0.075   0.9406    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.002793 on 192 degrees of freedom
## Multiple R-squared:  0.2179, Adjusted R-squared:  0.1975 
## F-statistic:  10.7 on 5 and 192 DF,  p-value: 4.392e-09

Visualization: Execution Stability

# Plot 1: Learning trend over time (Blocks 1-3)
p_exec_learning <- ggplot(exec_data_learning, aes(x = Total_Trial, y = Lyap_Acc, color = Sequence)) +
  geom_line(aes(group = interaction(Subject, Sequence)), alpha = 0.1) +
  # Use y ~ log(x) for logarithmic x-axis fitting
  geom_smooth(method = "lm", se = TRUE, formula = y ~ log(x)) + 
  #scale_x_log10(breaks = scales::trans_breaks("log10", function(x) 10^x),
  #              labels = scales::trans_format("log10", scales::math_format(10^.x))) + 
  labs(
    title = "Execution Stability Over Learning Phase (Blocks 1-3)",
    x = "Total Trial Number",
    y = expression("CoM Acceleration Lyapunov Exponent" ~ (Lyap["Acc"]))
  ) +
  theme_minimal(base_size = 14)

print(p_exec_learning)

# Plot 2: Boxplot comparison for Learning and Testing phases
p_exec_data <- ggplot(exec_data, aes(x = Phase, y = Lyap_Acc, fill = Sequence)) +
  geom_boxplot(alpha = 0.7) +
  labs(
    title = "Execution Stability Comparison (Learning vs Testing)",
    x = "Experimental Phase",
    y = expression("Lyapunov Exponent" ~ (Lyap["Acc"])),
    fill = "Sequence"
  ) +
  theme_minimal(base_size = 14)

print(p_exec_data)

Execution Variability vs RT

message("Linear Model for Execution Variability (RMS) vs. RT - Single Subject")
## Linear Model for Execution Variability (RMS) vs. RT - Single Subject
# LMEM: RMS_Acc ~ trial.RTS * Sequence + (1 | Subject)
# This model tests if trials with higher RT also have higher variability.

# NOTE: Since only one participant (PID) is loaded, we switch from LMEM (lmer) 
# to a simple Linear Model (lm) by removing the random effects term.
model_h5 <- lm(
  RMS_acc ~ trial.RTS * Sequence,
  data = exec_data
)

message("--- LM Summary for Execution RMS vs. RT ---")
## --- LM Summary for Execution RMS vs. RT ---
summary(model_h5)
## 
## Call:
## lm(formula = RMS_acc ~ trial.RTS * Sequence, data = exec_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.73078 -0.13600  0.03207  0.23930  0.94710 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            1.7254     0.1528  11.289  < 2e-16 ***
## trial.RTS             -0.4896     0.2089  -2.344  0.02010 *  
## Sequence12             0.7920     0.2134   3.711  0.00027 ***
## Sequence18             0.1395     0.2724   0.512  0.60902    
## trial.RTS:Sequence12  -0.7558     0.3517  -2.149  0.03289 *  
## trial.RTS:Sequence18   0.4375     0.4586   0.954  0.34133    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4557 on 192 degrees of freedom
## Multiple R-squared:  0.2834, Adjusted R-squared:  0.2647 
## F-statistic: 15.19 on 5 and 192 DF,  p-value: 1.449e-12