R Notebook

rm(list=ls(all=TRUE))  #clear previous

# Install pacman if not already installed
if (!require("pacman")) install.packages("pacman")

# Use pacman to load/install packages
pacman::p_load(
  gratia,
  gamm4,
  mgcv,
  segmented,
  nlme,
  mmrm,
  tidyverse,
  car,
  languageR,
  ggplot2,
  lattice,
  jtools,
  sjPlot,
  sjmisc,
  sjlabelled,
  gtools,
  lmerTest,
  ggeffects,
  ggpubr,
  corrplot,
  gridExtra,
  GGally,
  irr,
  blandr,
  knitr,
  kableExtra,
  ggprism,
  webshot2,
  htmltools,
  gt,
  parameters
)

# Read in full data, both pre and post assoc
AssocLongint <- read.csv(file = "/Volumes/vosslabhpc/Projects/BikeExtend/3-Experiment/2-Data/BIDS/derivatives/Bryan/Projects/EXTEND_beh/BikeEXTEND_Assoc_post/AssocPostSummary.csv")
AssocLongint$Subject <- as.factor(AssocLongint$Subject)
AssocLongint$Block <- as.numeric(AssocLongint$Block)
AssocLongint$Time <- as.factor(AssocLongint$Time)
AssocLongint$Intgrp <- ifelse(AssocLongint$Intgrp == "Group 1: Moderate Intensity", "ModVig", "Light"
)
AssocLongint <- AssocLongint[!is.na(AssocLongint$Intgrp), ]
AssocLongint$Intgrp <- as.factor(AssocLongint$Intgrp)


# Filter out pre data for Cross-sectional models
Assoc_pre <- AssocLongint %>%
  filter(Time == "preInt")

# Run GAMM model with only Block predicting Accuracy
gamm_model <- gamm(
  avg_acc ~ s(Block, k = 6),        # smooth term for block (nonlinear learning)
  random = list(Subject = ~1),      # random intercepts for each subject, 
  data = Assoc_pre
)

# ~1 indicates that each subject can have their own baseline (intercept), accounting for individual differences in average accuracy.

# This helps account for within-subject correlation, which is important in repeated-measures or longitudinal data (like multiple Blocks per subject)

# 'gamm()' fits a GAMM: it includes both smooth (nonlinear) terms and random effects.

# Unlike 'gam()', which only handles fixed effects, 'gamm()' is appropriate for repeated-measures data.

summary(gamm_model$gam)  # fixed effects and smooth term

Family: gaussian 
Link function: identity 

Formula:
avg_acc ~ s(Block, k = 6)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.47901    0.01464   32.72   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
           edf Ref.df     F p-value    
s(Block) 3.262  3.262 486.7  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.294   
  Scale est. = 0.010265  n = 1150

summary(gamm_model$lme)  # random effects (subject-level)

Linear mixed-effects model fit by maximum likelihood
  Data: strip.offset(mf) 

Random effects:
 Formula: ~Xr - 1 | g
 Structure: pdIdnot
               Xr1        Xr2        Xr3        Xr4
StdDev: 0.06800683 0.06800683 0.06800683 0.06800683

 Formula: ~1 | Subject %in% g
        (Intercept)  Residual
StdDev:   0.1536245 0.1013163

Fixed effects:  y ~ X - 1 

 Correlation: 
             X(Int)
Xs(Block)Fx1 0     

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-3.02328817 -0.60989638 -0.01281136  0.64526397  3.80618117 

Number of Observations: 1150
Number of Groups: 
             g Subject %in% g 
             1            115 

plot(gamm_model$gam, residuals = TRUE, pch = 16, main = "Smooth Term: ACC ~ BlockNum")

#library(gratia)
draw(gamm_model$gam) # uses ggplot under the hood for more control over aesthetics

# gamm4() uses the lme4 engine for random effects, which is more flexible and scalable than mgcv::gamm()


gamm4_model <- gamm4(
  avg_acc ~ s(Block, k = 6),     # Fixed effect: smooth spline over Block to model nonlinear learning
  random = ~(1 | Subject),       # Random effect: random intercepts for each subject
  data = Assoc_pre              
)

# Summary of the GAM component: shows fixed effects and smooth terms
summary(gamm4_model$gam)

Family: gaussian 
Link function: identity 

Formula:
avg_acc ~ s(Block, k = 6)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   0.4790     0.0147   32.59   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
           edf Ref.df     F p-value    
s(Block) 3.503  3.503 453.7  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.294   
lmer.REML = -1615.5  Scale est. = 0.010265  n = 1150

# Summary of the Mixed-Effects component: shows variance due to subjects (random effects)
summary(gamm4_model$lme)

Length  Class   Mode 
     0   NULL   NULL 

# Diagnostic plot: smooth function with residuals overlaid
plot(gamm4_model$gam, residuals = TRUE, pch = 16)

# Fit a Generalized Additive Mixed Model (GAMM) using mgcv::gamm()

gamm_model_group <- gamm(
  avg_acc ~ Intgrp + Time + Age + Sex + YrsEdu +           # Parametric fixed effects
    s(Block, by = interaction(Intgrp, Time), k = 6),       # Group × Time-specific nonlinear learning curve
  random = list(Subject = ~1),                             # Random intercepts for each subject (nested design)
  data = AssocLongint                                      # Long-format dataset including pre and post data
)

# View summary of the GAM portion (fixed effects + smooth terms)
summary(gamm_model_group$gam)

Family: gaussian 
Link function: identity 

Formula:
avg_acc ~ Intgrp + Time + Age + Sex + YrsEdu + s(Block, by = interaction(Intgrp, 
    Time), k = 6)

Parametric coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)   1.015223   0.151990   6.680 3.05e-11 ***
IntgrpModVig -0.059552   0.025638  -2.323   0.0203 *  
TimepreInt   -0.005090   0.005310  -0.959   0.3379    
Age          -0.010134   0.002286  -4.433 9.77e-06 ***
SexMale      -0.012511   0.027061  -0.462   0.6439    
YrsEdu        0.008795   0.005598   1.571   0.1163    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
                                                   edf Ref.df     F p-value    
s(Block):interaction(Intgrp, Time)Light.postInt  2.347  2.347 264.9  <2e-16 ***
s(Block):interaction(Intgrp, Time)ModVig.postInt 3.776  3.776 120.8  <2e-16 ***
s(Block):interaction(Intgrp, Time)Light.preInt   2.971  2.971 197.6  <2e-16 ***
s(Block):interaction(Intgrp, Time)ModVig.preInt  2.082  2.082 281.0  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.367   
  Scale est. = 0.013856  n = 2140

# View summary of the LME portion (random effects for subjects)
summary(gamm_model_group$lme)

Linear mixed-effects model fit by maximum likelihood
  Data: strip.offset(mf) 

Random effects:
 Formula: ~Xr - 1 | g
 Structure: pdIdnot
               Xr1        Xr2        Xr3        Xr4
StdDev: 0.05642862 0.05642862 0.05642862 0.05642862

 Formula: ~Xr.0 - 1 | g.0 %in% g
 Structure: pdIdnot
            Xr.01     Xr.02     Xr.03     Xr.04
StdDev: 0.1797177 0.1797177 0.1797177 0.1797177

 Formula: ~Xr.1 - 1 | g.1 %in% g.0 %in% g
 Structure: pdIdnot
             Xr.11      Xr.12      Xr.13      Xr.14
StdDev: 0.09001205 0.09001205 0.09001205 0.09001205

 Formula: ~Xr.2 - 1 | g.2 %in% g.1 %in% g.0 %in% g
 Structure: pdIdnot
             Xr.21      Xr.22      Xr.23      Xr.24
StdDev: 0.03975747 0.03975747 0.03975747 0.03975747

 Formula: ~1 | Subject %in% g.2 %in% g.1 %in% g.0 %in% g
        (Intercept)  Residual
StdDev:   0.1345361 0.1177126

Fixed effects:  y ~ X - 1 

 Correlation: 
                                                     X(Int) XIntMV XTmprI XAge   XSexMl XYrsEd Xs(Blck):ntrctn(Intgrp,Tm)Lght.psIF1 Xs(Blck):ntrctn(Intgrp,Tm)MdVg.psIF1 Xs(Blck):ntrctn(Intgrp,Tm)Lght.prIF1
XIntgrpModVig                                        -0.108                                                                                                                                                  
XTimepreInt                                          -0.025 -0.003                                                                                                                                           
XAge                                                 -0.798  0.000  0.005                                                                                                                                    
XSexMale                                              0.006  0.059 -0.010 -0.029                                                                                                                             
XYrsEdu                                              -0.368  0.029  0.001 -0.251 -0.074                                                                                                                      
Xs(Block):interaction(Intgrp, Time)Light.postIntFx1   0.000  0.000  0.000  0.000  0.000  0.000                                                                                                               
Xs(Block):interaction(Intgrp, Time)ModVig.postIntFx1  0.000  0.000  0.000  0.000  0.000  0.000  0.000                                                                                                        
Xs(Block):interaction(Intgrp, Time)Light.preIntFx1    0.000  0.000  0.000  0.000  0.000  0.000  0.000                                0.000                                                                   
Xs(Block):interaction(Intgrp, Time)ModVig.preIntFx1   0.000  0.000  0.000  0.000  0.000  0.000  0.000                                0.000                                0.000                              

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-3.55827765 -0.63617325  0.03270191  0.68153449  3.66596309 

Number of Observations: 2140
Number of Groups: 
                                        g                                g.0 %in% g                       g.1 %in% g.0 %in% g              g.2 %in% g.1 %in% g.0 %in% g Subject %in% g.2 %in% g.1 %in% g.0 %in% g 
                                        1                                         1                                         1                                         1                                       116 

# Plot the smooth learning curves for each Group × Time interaction
# Shows the estimated smooth function with 95% confidence bands and residuals
plot(gamm_model_group$gam, pages = 1, residuals = TRUE, pch = 16)

NA
NA

# Fit a Generalized Additive Mixed Model using gamm4 (which uses lme4::lmer for random effects). Repeated model from above that used gamm()
gamm4_model_group <- gamm4(
  avg_acc ~ Intgrp + Time + Age + Sex + YrsEdu +          # Parametric fixed effects
    s(Block, by = interaction(Intgrp, Time), k = 6),      # Smooth term for Block by Group × Time interaction:Fits a separate smooth learning curve over Block for each unique combination of intervention group (Intgrp) and                                                                   timepoint (Time).The interaction() function defines these subgroups.
  random = ~(1 | Subject),                                # Random intercepts for each subject (lme4-style formula)
  data = AssocLongint                                     
)

# Summarize fixed effects and smooth terms from the GAM part (mgcv::gam object)
summary(gamm4_model_group$gam)

Family: gaussian 
Link function: identity 

Formula:
avg_acc ~ Intgrp + Time + Age + Sex + YrsEdu + s(Block, by = interaction(Intgrp, 
    Time), k = 6)

Parametric coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)   1.015231   0.155175   6.542 7.56e-11 ***
IntgrpModVig -0.059574   0.026175  -2.276   0.0229 *  
TimepreInt   -0.005087   0.005305  -0.959   0.3378    
Age          -0.010135   0.002334  -4.342 1.48e-05 ***
SexMale      -0.012517   0.027627  -0.453   0.6506    
YrsEdu        0.008796   0.005715   1.539   0.1239    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
                                                   edf Ref.df     F p-value    
s(Block):interaction(Intgrp, Time)Light.postInt  2.587  2.587 241.0  <2e-16 ***
s(Block):interaction(Intgrp, Time)ModVig.postInt 3.896  3.896 117.7  <2e-16 ***
s(Block):interaction(Intgrp, Time)Light.preInt   3.259  3.259 180.8  <2e-16 ***
s(Block):interaction(Intgrp, Time)ModVig.preInt  2.403  2.403 244.3  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.367   
lmer.REML = -2618.3  Scale est. = 0.013867  n = 2140

# Summarize random effects from the LME part (lme4::merMod object)
summary(gamm4_model_group$mer)

Linear mixed model fit by REML ['lmerMod']

REML criterion at convergence: -2618.3

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.5568 -0.6334  0.0306  0.6855  3.6840 

Random effects:
 Groups   Name                                             Variance Std.Dev.
 Subject  (Intercept)                                      0.018954 0.13767 
 Xr.2     s(Block):interaction(Intgrp, Time)ModVig.preInt  0.002937 0.05419 
 Xr.1     s(Block):interaction(Intgrp, Time)Light.preInt   0.012783 0.11306 
 Xr.0     s(Block):interaction(Intgrp, Time)ModVig.postInt 0.039432 0.19858 
 Xr       s(Block):interaction(Intgrp, Time)Light.postInt  0.004896 0.06997 
 Residual                                                  0.013867 0.11776 
Number of obs: 2140, groups:  Subject, 116; Xr.2, 4; Xr.1, 4; Xr.0, 4; Xr, 4

Fixed effects:
                                                      Estimate Std. Error t value
X(Intercept)                                          1.015231   0.155175   6.542
XIntgrpModVig                                        -0.059574   0.026175  -2.276
XTimepreInt                                          -0.005087   0.005305  -0.959
XAge                                                 -0.010135   0.002334  -4.342
XSexMale                                             -0.012517   0.027627  -0.453
XYrsEdu                                               0.008796   0.005715   1.539
Xs(Block):interaction(Intgrp, Time)Light.postIntFx1  -0.139448   0.020447  -6.820
Xs(Block):interaction(Intgrp, Time)ModVig.postIntFx1 -0.196712   0.034207  -5.751
Xs(Block):interaction(Intgrp, Time)Light.preIntFx1   -0.147200   0.025848  -5.695
Xs(Block):interaction(Intgrp, Time)ModVig.preIntFx1  -0.119246   0.016805  -7.096

Correlation of Fixed Effects:
                                     X(Int) XIntMV XTmprI XAge   XSexMl XYrsEd Xs(Blck):ntrctn(Intgrp,Tm)Lght.psIF1 Xs(Blck):ntrctn(Intgrp,Tm)MdVg.psIF1 Xs(Blck):ntrctn(Intgrp,Tm)Lght.prIF1
XIntgrpMdVg                          -0.108                                                                                                                                                  
XTimepreInt                          -0.024 -0.003                                                                                                                                           
XAge                                 -0.798  0.000  0.005                                                                                                                                    
XSexMale                              0.006  0.059 -0.009 -0.029                                                                                                                             
XYrsEdu                              -0.368  0.029  0.001 -0.251 -0.074                                                                                                                      
Xs(Blck):ntrctn(Intgrp,Tm)Lght.psIF1  0.000  0.000  0.000  0.000  0.000  0.000                                                                                                               
Xs(Blck):ntrctn(Intgrp,Tm)MdVg.psIF1  0.000  0.000  0.000  0.000  0.000  0.000  0.000                                                                                                        
Xs(Blck):ntrctn(Intgrp,Tm)Lght.prIF1  0.000  0.000  0.000  0.000  0.000  0.000  0.000                                0.000                                                                   
Xs(Blck):ntrctn(Intgrp,Tm)MdVg.prIF1  0.000  0.000  0.000  0.000  0.000  0.000  0.000                                0.000                                0.000                              

# Plot the smooth terms with residuals
plot(gamm4_model_group$gam, residuals = TRUE, pch = 16)

# Step 1: Extract smooth term estimates from the fitted GAMM model
# This includes estimated values and standard errors for each smooth function in the model
smooths_df <- gratia::smooth_estimates(gamm_model_group$gam)

# Step 2: Inspect the column names to identify how smooth terms are labeled
# Useful for confirming which column contains the smooth labels (e.g., ".smooth")
print(colnames(smooths_df))  # Diagnostic step to understand structure of smooths_df

[1] ".smooth"                   ".type"                     ".by"                       ".estimate"                 ".se"                       "Block"                     "interaction(Intgrp, Time)"

# Step 3: Create label for each group × time interaction
# Based on pattern matching within the `.smooth` column
# This is used for coloring and grouping in the plot
smooths_df <- smooths_df %>%
  mutate(Group = case_when(
    grepl("Light\\.postInt", .smooth) ~ "Light - Post",
    grepl("ModVig\\.postInt", .smooth) ~ "ModVig - Post",
    grepl("Light\\.preInt", .smooth) ~ "Light - Pre",
    grepl("ModVig\\.preInt", .smooth) ~ "ModVig - Pre",
    TRUE ~ "Other"  # Catch-all fallback (shouldn’t typically occur)
  ))

# Step 1: Extract smooth term estimates from the fitted GAMM model

There were 50 or more warnings (use warnings() to see the first 50)

# This includes estimated values and standard errors for each smooth function in the model
smooths_df4 <- gratia::smooth_estimates(gamm4_model_group$gam)

# Step 2: Inspect the column names to identify how smooth terms are labeled
# Useful for confirming which column contains the smooth labels (e.g., ".smooth")
print(colnames(smooths_df4))  # Diagnostic step to understand structure of smooths_df

[1] ".smooth"                   ".type"                     ".by"                       ".estimate"                 ".se"                       "Block"                     "interaction(Intgrp, Time)"

# Step 3: Create label for each group × time interaction
# Based on pattern matching within the `.smooth` column
# This is used for coloring and grouping in the plot
smooths_df4 <- smooths_df4 %>%
  mutate(Group = case_when(
    grepl("Light\\.postInt", .smooth) ~ "Light - Post",
    grepl("ModVig\\.postInt", .smooth) ~ "ModVig - Post",
    grepl("Light\\.preInt", .smooth) ~ "Light - Pre",
    grepl("ModVig\\.preInt", .smooth) ~ "ModVig - Pre",
    TRUE ~ "Other"  # Catch-all fallback (shouldn’t typically occur)
  ))

# Step 4: Create a TVEM-style line plot showing estimated learning curves over Block
# Each line represents the smooth trajectory for a Group × Time combination
# A ribbon is added to show ±1 standard error around the estimated smooth
ggplot(smooths_df, aes(x = Block, y = .estimate, color = Group, fill = Group)) +
  geom_line(size = 1.2) +  # Line plot of estimated smooth
  geom_ribbon(aes(ymin = .estimate - .se, ymax = .estimate + .se), 
              alpha = 0.2, color = NA) +  # Shaded confidence band
  labs(
    title = "TVEM-style Plot of Accuracy by Block and Group × Time\n - Gamm",
    x = "Block", y = "Estimated Accuracy",
    color = "Group × Time", fill = "Group × Time"
  ) +
  theme_minimal(base_size = 14)  # Clean visual style

ggplot(smooths_df4, aes(x = Block, y = .estimate, color = Group, fill = Group)) +
  geom_line(size = 1.2) +  # Line plot of estimated smooth
  geom_ribbon(aes(ymin = .estimate - .se, ymax = .estimate + .se), 
              alpha = 0.2, color = NA) +  # Shaded confidence band
  labs(
    title = "TVEM-style Plot of Accuracy by Block and Group × Time\n - Gamm4",
    x = "Block", y = "Estimated Accuracy",
    color = "Group × Time", fill = "Group × Time"
  ) +
  theme_minimal(base_size = 14)  # Clean visual style

# Step 1: Add TimePeriod column
smooths_df <- smooths_df %>%
  mutate(TimePeriod = case_when(
    grepl("preInt", .smooth) ~ "Pre",
    grepl("postInt", .smooth) ~ "Post",
    TRUE ~ "Other"
  ))

# Step 2: Plot with facets
ggplot(smooths_df, aes(x = Block, y = .estimate, color = Group, fill = Group)) +
  geom_line(size = 1.2) +
  geom_ribbon(aes(ymin = .estimate - .se, ymax = .estimate + .se), alpha = 0.2, color = NA) +
  facet_wrap(~TimePeriod) +
  labs(
    title = "TVEM-style Plot of Accuracy by Block - Gamm",
    x = "Block", y = "Estimated Accuracy",
    color = "Group", fill = "Group"
  ) +
  theme_minimal(base_size = 14)

# Step 1: Add TimePeriod column
smooths_df4 <- smooths_df4 %>%
  mutate(TimePeriod = case_when(
    grepl("preInt", .smooth) ~ "Pre",
    grepl("postInt", .smooth) ~ "Post",
    TRUE ~ "Other"
  ))

# Step 2: Plot with facets
ggplot(smooths_df4, aes(x = Block, y = .estimate, color = Group, fill = Group)) +
  geom_line(size = 1.2) +
  geom_ribbon(aes(ymin = .estimate - .se, ymax = .estimate + .se), alpha = 0.2, color = NA) +
  facet_wrap(~TimePeriod) +
  labs(
    title = "TVEM-style Plot of Accuracy by Block - Gamm4",
    x = "Block", y = "Estimated Accuracy",
    color = "Group", fill = "Group"
  ) +
  theme_minimal(base_size = 14)

#

#

`

This section look at the EXTEND baseline data with the addition of the Midlife sample

These demographics has only Age for midlife because at the time collected, rest of demo var were missing. Need to return to Redcap to retrieve the rest

# TaskSwitchMOb <- read.csv("/Volumes/vosslabhpc/Projects/BikeExtend/3-Experiment/2-Data/Participants/BaselineSummaries/TaskSwitch_Trial_Pre_Full.csv", stringsAsFactors = FALSE)[, -1]
# TaskSwitchMOb <- TaskSwitchMOb[, c(1:7,12:33)]

# This df has already dropped the same subs dropped in the intervention analysis
TaskSwitchMO <- read.csv("/Volumes/vosslabhpc/Projects/BikeExtend/3-Experiment/2-Data/BIDS/derivatives/Bryan/Projects/Task_Switch-Intervention/Arch/Pre_TS_Long.csv", stringsAsFactors = FALSE)[, -1] %>% mutate(Time = "pre")
TaskSwitchMO <- TaskSwitchMO[, -c(31:35)]
bike <- read.csv("/Volumes/vosslabhpc/Projects/BikeExtend/3-Experiment/2-Data/BIDS/derivatives/Bryan/Projects/TertiarySulci/EXTEND_labels/code/00_Analysis/00_Variables/BikeDemo.csv", header = TRUE, sep = ",", stringsAsFactors = FALSE)
colnames(bike)[1:12] <- c("Subject", "Sex", "Age", "YrsEdu", "Moca", "RAVLT", "RAVLTpost",
                          "AnimalFlency", "Intgrp", "Rel_Vo2m", "RespExRatio", "DSST")
bike <- bike[, -c(2,3)]

# Drop rows with missing cardiorespiratory fitness data
bike <- bike %>% drop_na(Rel_Vo2m)

TaskSwitchMO <- TaskSwitchMO %>%
  left_join(bike, by = "Subject")

# Ensure Subject is treated as a character to extract prefix to create grouping variable
TaskSwitchMO <- TaskSwitchMO %>%
  mutate(
    Subject = as.character(Subject),
    EXTGroup = case_when(
      str_starts(Subject, "5") ~ "Midlife",
      str_starts(Subject, "2") ~ "Older",
      TRUE ~ NA_character_  # in case there are other IDs
    )
  )

TaskSwitchMO$Subject <- as.factor(TaskSwitchMO$Subject)
TaskSwitchMO$EXTGroup <- as.factor(TaskSwitchMO$EXTGroup)
TaskSwitchMO$Age <- as.integer(TaskSwitchMO$Age)
TaskSwitchMO$YrsEdu <- as.integer(TaskSwitchMO$YrsEdu)
TaskSwitchMO$Condition <- as.factor(TaskSwitchMO$Condition)
TaskSwitchMO$Trial.Type <- as.factor(TaskSwitchMO$Trial.Type)
TaskSwitchMO$Sex <- as.factor(TaskSwitchMO$Sex)

contrasts(TaskSwitchMO$Trial.Type)

         Single Switched
Repeated      0        0
Single        1        0
Switched      0        1

print("Mixing Cost is referencing Repaeted, Single-Repeat, therefore take the negative of the beta shown for single (multiply by -1)")

[1] "Mixing Cost is referencing Repaeted, Single-Repeat, therefore take the negative of the beta shown for single (multiply by -1)"

print("Switching Cost is referenced here correctly, Switching - Repeated, therfore interpret the beta as is")

[1] "Switching Cost is referenced here correctly, Switching - Repeated, therfore interpret the beta as is"

Analysing Age as a continous variable

Age + Condition * Trial.Type

Age <- Min=40, Max=80

Agingsimple_model <- lmer(
  Correct.RT ~ Age + Condition * Trial.Type + (1 | Subject),
  data = TaskSwitchMO
)
summary(Agingsimple_model) #, ddf=c("Kenward-Roger"))

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Correct.RT ~ Age + Condition * Trial.Type + (1 | Subject)
   Data: TaskSwitchMO

REML criterion at convergence: 1215709

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.0743 -0.6255 -0.1863  0.4148  5.5438 

Random effects:
 Groups   Name        Variance Std.Dev.
 Subject  (Intercept)  12023   109.7   
 Residual             110000   331.7   
Number of obs: 84115, groups:  Subject, 148

Fixed effects:
                                      Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)                            913.114     76.403   156.049  11.951  < 2e-16 ***
Age                                      3.596      1.237   155.862   2.908  0.00417 ** 
ConditionSeparate                     -373.395      3.892 83999.363 -95.936  < 2e-16 ***
Trial.TypeSingle                      -287.071      4.042 83966.690 -71.020  < 2e-16 ***
Trial.TypeSwitched                     232.757      4.049 83967.964  57.483  < 2e-16 ***
ConditionSeparate:Trial.TypeSingle     275.732      5.608 83964.760  49.165  < 2e-16 ***
ConditionSeparate:Trial.TypeSwitched   -57.726      5.559 83966.107 -10.384  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
                    (Intr) Age    CndtnS Trl.TypSn Trl.TypSw CndtnSprt:Trl.TypSn
Age                 -0.992                                                      
ConditnSprt         -0.026 -0.001                                               
Trl.TypSngl         -0.025 -0.001  0.501                                        
Trl.TypSwtc         -0.027  0.001  0.500  0.482                                 
CndtnSprt:Trl.TypSn  0.018  0.000 -0.689 -0.720    -0.347                       
CndtnSprt:Trl.TypSw  0.020 -0.001 -0.694 -0.351    -0.727     0.482             

print("Note: Recall to interpret the single trial betas as negative (multiply by -1) to infer Mixing Cost")

[1] "Note: Recall to interpret the single trial betas as negative (multiply by -1) to infer Mixing Cost"

# Predict RT across Condition and Trial.Type for three representative age values
Agingsimple_model.predict <- ggpredict(
  Agingsimple_model,
  terms = c("Condition", "Trial.Type")  # You can change these age values as needed
)

Agingsimple_model.PLOT <- ggplot(Agingsimple_model.predict, aes(x = x, y = predicted, group = group)) +
  geom_line(aes(linetype = group, color = group), size = 1.2) +  # Trial.Type as line types/colors
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = group), 
              alpha = 0.2, color = NA) +
  theme_bw() +
  theme(
    text = element_text(size = 16),
    strip.background = element_blank(),
    strip.text = element_text(size = 12, face = "bold"),
    legend.position = "bottom",
    legend.title = element_text(size = 11),
    legend.text = element_text(size = 11),
    plot.title = element_text(hjust = 0.5)
  ) +
  labs(
    title = "Predicted RT by Condition and Trial Type",
    x = "Condition",
    y = "Predicted RT (ms)",
    linetype = "Trial Type",
    color = "Trial Type",
    fill = "Trial Type"
  )
Agingsimple_model.PLOT

Analysing Age as a continous variable

Age + Condition * Trial.Type

Age <- Min=40, Max=80

Agingsimple_model <- lmer(
  Correct.RT ~ Age + Condition * Trial.Type + (1 | Subject),
  data = TaskSwitchMO
)
summary(Agingsimple_model) #, ddf=c("Kenward-Roger"))

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Correct.RT ~ Age + Condition * Trial.Type + (1 | Subject)
   Data: TaskSwitchMO

REML criterion at convergence: 1215709

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.0743 -0.6255 -0.1863  0.4148  5.5438 

Random effects:
 Groups   Name        Variance Std.Dev.
 Subject  (Intercept)  12023   109.7   
 Residual             110000   331.7   
Number of obs: 84115, groups:  Subject, 148

Fixed effects:
                                      Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)                            913.114     76.403   156.049  11.951  < 2e-16 ***
Age                                      3.596      1.237   155.862   2.908  0.00417 ** 
ConditionSeparate                     -373.395      3.892 83999.363 -95.936  < 2e-16 ***
Trial.TypeSingle                      -287.071      4.042 83966.690 -71.020  < 2e-16 ***
Trial.TypeSwitched                     232.757      4.049 83967.964  57.483  < 2e-16 ***
ConditionSeparate:Trial.TypeSingle     275.732      5.608 83964.760  49.165  < 2e-16 ***
ConditionSeparate:Trial.TypeSwitched   -57.726      5.559 83966.107 -10.384  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
                    (Intr) Age    CndtnS Trl.TypSn Trl.TypSw CndtnSprt:Trl.TypSn
Age                 -0.992                                                      
ConditnSprt         -0.026 -0.001                                               
Trl.TypSngl         -0.025 -0.001  0.501                                        
Trl.TypSwtc         -0.027  0.001  0.500  0.482                                 
CndtnSprt:Trl.TypSn  0.018  0.000 -0.689 -0.720    -0.347                       
CndtnSprt:Trl.TypSw  0.020 -0.001 -0.694 -0.351    -0.727     0.482             

print("Note: Recall to interpret the single trial betas as negative (multiply by -1) to infer Mixing Cost")

[1] "Note: Recall to interpret the single trial betas as negative (multiply by -1) to infer Mixing Cost"

# Predict RT across Condition and Trial.Type for three representative age values
Agingsimple_model.predict <- ggpredict(
  Agingsimple_model,
  terms = c("Condition", "Trial.Type")  # You can change these age values as needed
)

Agingsimple_model.PLOT <- ggplot(Agingsimple_model.predict, aes(x = x, y = predicted, group = group)) +
  geom_line(aes(linetype = group, color = group), size = 1.2) +  # Trial.Type as line types/colors
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = group), 
              alpha = 0.2, color = NA) +
  theme_bw() +
  theme(
    text = element_text(size = 16),
    strip.background = element_blank(),
    strip.text = element_text(size = 12, face = "bold"),
    legend.position = "bottom",
    legend.title = element_text(size = 11),
    legend.text = element_text(size = 11),
    plot.title = element_text(hjust = 0.5)
  ) +
  labs(
    title = "Predicted RT by Condition and Trial Type at Selected Ages",
    x = "Condition",
    y = "Predicted RT (ms)",
    linetype = "Trial Type",
    color = "Trial Type",
    fill = "Trial Type"
  )
Agingsimple_model.PLOT

# Predict RT across Condition and Trial.Type for three representative age values
Agingsimple_model.predict <- ggpredict(
  Agingsimple_model,
  terms = c("Condition", "Trial.Type", "Age [50,65,80]")  # You can change these age values as needed
)

Agingsimple.PLOT2 <- ggplot(Agingsimple_model.predict, aes(x = x, y = predicted, group = group)) +
  geom_line(aes(linetype = group, color = group), size = 1.2) +  # Trial.Type as line types/colors
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = group), 
              alpha = 0.2, color = NA) +
  facet_wrap(~facet, labeller = label_both) +  # Facets by Age value
  ylim(600, 1500) +
  theme_bw() +
  theme(
    text = element_text(size = 16),
    strip.background = element_blank(),
    strip.text = element_text(size = 12, face = "bold"),
    legend.position = "bottom",
    legend.title = element_text(size = 11),
    legend.text = element_text(size = 11),
    plot.title = element_text(hjust = 0.5)
  ) +
  labs(
    title = "Predicted RT by Condition and Trial Type\n at Selected Ages",
    x = "Condition",
    y = "Predicted RT (ms)",
    linetype = "Trial Type",
    color = "Trial Type",
    fill = "Trial Type"
  )
Agingsimple.PLOT2

Analysing Age as a continous variable

Age + Condition * Trial.Type

Age <- Min=40, Max=80

Fitness_simple_model <- lmer(
  Correct.RT ~ Age + Rel_Vo2m + Condition * Trial.Type + (1 | Subject),
  data = TaskSwitchMO
)
summary(Fitness_simple_model) #, ddf=c("Kenward-Roger"))

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Correct.RT ~ Age + Rel_Vo2m + Condition * Trial.Type + (1 | Subject)
   Data: TaskSwitchMO

REML criterion at convergence: 961599.8

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.0723 -0.6276 -0.1876  0.4193  5.4583 

Random effects:
 Groups   Name        Variance Std.Dev.
 Subject  (Intercept)  11829   108.8   
 Residual             110720   332.7   
Number of obs: 66504, groups:  Subject, 118

Fixed effects:
                                      Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)                            957.010    130.145   115.148   7.353 3.05e-11 ***
Age                                      3.789      1.788   115.080   2.119   0.0362 *  
Rel_Vo2m                                -2.817      1.957   114.859  -1.439   0.1528    
ConditionSeparate                     -368.451      4.386 66386.032 -83.999  < 2e-16 ***
Trial.TypeSingle                      -280.990      4.583 66385.389 -61.312  < 2e-16 ***
Trial.TypeSwitched                     241.920      4.602 66386.341  52.564  < 2e-16 ***
ConditionSeparate:Trial.TypeSingle     268.785      6.332 66383.294  42.450  < 2e-16 ***
ConditionSeparate:Trial.TypeSwitched   -63.766      6.282 66384.363 -10.151  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
                    (Intr) Age    Rl_V2m CndtnS Trl.TypSn Trl.TypSw CndtnSprt:Trl.TypSn
Age                 -0.950                                                             
Rel_Vo2m            -0.536  0.258                                                      
ConditnSprt         -0.017 -0.001  0.001                                               
Trl.TypSngl         -0.017  0.000  0.000  0.506                                        
Trl.TypSwtc         -0.018  0.001 -0.001  0.503  0.483                                 
CndtnSprt:Trl.TypSn  0.012  0.000  0.000 -0.692 -0.723    -0.349                       
CndtnSprt:Trl.TypSw  0.013 -0.001  0.000 -0.697 -0.353    -0.732     0.483             

print("Note: Recall to interpret the single trial betas as negative (multiply by -1) to infer Mixing Cost")

[1] "Note: Recall to interpret the single trial betas as negative (multiply by -1) to infer Mixing Cost"

# Predict RT across Condition and Trial.Type for three representative age values
Fitness_simple_model.predict <- ggpredict(
  Fitness_simple_model,
  terms = c("Condition", "Trial.Type")  # You can change these age values as needed
)

Fitness_simple_model.PLOT <- ggplot(Fitness_simple_model.predict, aes(x = x, y = predicted, group = group)) +
  geom_line(aes(linetype = group, color = group), size = 1.2) +  # Trial.Type as line types/colors
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = group), 
              alpha = 0.2, color = NA) +
  theme_bw() +
  theme(
    text = element_text(size = 16),
    strip.background = element_blank(),
    strip.text = element_text(size = 12, face = "bold"),
    legend.position = "bottom",
    legend.title = element_text(size = 11),
    legend.text = element_text(size = 11),
    plot.title = element_text(hjust = 0.5)
  ) +
  labs(
    title = "Predicted RT by Condition and Trial Type",
    x = "Condition",
    y = "Predicted RT (ms)",
    linetype = "Trial Type",
    color = "Trial Type",
    fill = "Trial Type"
  )
Fitness_simple_model.PLOT

# Predict RT across Condition and Trial.Type for three representative age values
Fitness_simple_model.predict <- ggpredict(
  Fitness_simple_model,
  terms = c("Condition", "Trial.Type", "Rel_Vo2m [16.30,19.70,24.90]")  # You can change these age values as needed
)

Fitness_simple_model.PLOT2 <- ggplot(Fitness_simple_model.predict, aes(x = x, y = predicted, group = group)) +
  geom_line(aes(linetype = group, color = group), size = 1.2) +  # Trial.Type as line types/colors
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = group), 
              alpha = 0.2, color = NA) +
  facet_wrap(~facet, labeller = label_both) +  # Facets by Age value
  ylim(600, 1500) +
  theme_bw() +
  theme(
    text = element_text(size = 16),
    strip.background = element_blank(),
    strip.text = element_text(size = 12, face = "bold"),
    legend.position = "bottom",
    legend.title = element_text(size = 11),
    legend.text = element_text(size = 11),
    plot.title = element_text(hjust = 0.5)
  ) +
  labs(
    title = "Predicted RT by Condition and Trial Type\n at Selected RelVo2m",
    x = "Condition",
    y = "Predicted RT (ms)",
    linetype = "Trial Type",
    color = "Trial Type",
    fill = "Trial Type"
  )
Fitness_simple_model.PLOT2

Fitness_simple_model.PLOT2

Agingsimple.PLOT2

Analysing Age as a continous variable

Age * Condition * Trial.Type

Age <- Min=40, Max=80

Analysing Age as a continous variable w/ the agegroup variable (Midlife and Older)

Condition * Trial.Type * AgeGroup - with Age as a covariate

Age <- Min=40, Max=80