MasterThesis_Analysis

Packages

library(tidyverse)

Warning: package 'tidyverse' was built under R version 4.2.3

Warning: package 'tidyr' was built under R version 4.2.3

Warning: package 'readr' was built under R version 4.2.3

Warning: package 'dplyr' was built under R version 4.2.3

Warning: package 'stringr' was built under R version 4.2.3

Warning: package 'forcats' was built under R version 4.2.3

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.2     ✔ 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(readxl)
library(tidyr)
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(effects)

Warning: package 'effects' was built under R version 4.2.3

Loading required package: carData
lattice theme set by effectsTheme()
See ?effectsTheme for details.

library(emmeans)

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

library(haven)
library(effects)
library(lattice)
library(car)

Attaching package: 'car'

The following object is masked from 'package:dplyr':

    recode

The following object is masked from 'package:purrr':

    some

library(ggplot2)
library(knitr)
library(reshape2)

Attaching package: 'reshape2'

The following object is masked from 'package:tidyr':

    smiths

library(dplyr)
library(forcats)
library(DHARMa)

This is DHARMa 0.4.7. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa')

library(Hmisc)

Attaching package: 'Hmisc'

The following objects are masked from 'package:dplyr':

    src, summarize

The following objects are masked from 'package:base':

    format.pval, units

library(phia)
library(lsmeans)

The 'lsmeans' package is now basically a front end for 'emmeans'.
Users are encouraged to switch the rest of the way.
See help('transition') for more information, including how to
convert old 'lsmeans' objects and scripts to work with 'emmeans'.

library(multcomp)

Loading required package: mvtnorm
Loading required package: survival
Loading required package: TH.data
Loading required package: MASS

Attaching package: 'MASS'

The following object is masked from 'package:dplyr':

    select


Attaching package: 'TH.data'

The following object is masked from 'package:MASS':

    geyser

library(plotly)

Attaching package: 'plotly'

The following object is masked from 'package:MASS':

    select

The following object is masked from 'package:Hmisc':

    subplot

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

library(lmerTest)
library(tinytex)
library(devtools)

Loading required package: usethis

Attaching package: 'devtools'

The following object is masked from 'package:emmeans':

    test

library(brms)

Loading required package: Rcpp
Loading 'brms' package (version 2.22.0). Useful instructions
can be found by typing help('brms'). A more detailed introduction
to the package is available through vignette('brms_overview').

Attaching package: 'brms'

The following object is masked from 'package:survival':

    kidney

The following object is masked from 'package:lme4':

    ngrps

The following object is masked from 'package:stats':

    ar

library(janitor)

Attaching package: 'janitor'

The following objects are masked from 'package:stats':

    chisq.test, fisher.test

library(magrittr)

Attaching package: 'magrittr'

The following object is masked from 'package:purrr':

    set_names

The following object is masked from 'package:tidyr':

    extract

library(mascutils)

Attaching package: 'mascutils'

The following object is masked from 'package:car':

    logit

The following object is masked from 'package:tidyr':

    expand_grid

The following object is masked from 'package:base':

    mode

library(ggthemes)

Warning: package 'ggthemes' was built under R version 4.2.3

library(afex)

************
Welcome to afex. For support visit: http://afex.singmann.science/
- Functions for ANOVAs: aov_car(), aov_ez(), and aov_4()
- Methods for calculating p-values with mixed(): 'S', 'KR', 'LRT', and 'PB'
- 'afex_aov' and 'mixed' objects can be passed to emmeans() for follow-up tests
- Get and set global package options with: afex_options()
- Set sum-to-zero contrasts globally: set_sum_contrasts()
- For example analyses see: browseVignettes("afex")
************

Attaching package: 'afex'

The following object is masked from 'package:lme4':

    lmer

library(openxlsx)

Importing data

demographics <- read_excel("PilotPerformance.xlsx", sheet = 2)
nasa <- read_excel("PilotPerformance.xlsx", sheet = 3)

New names:
• `` -> `...13`
• `` -> `...14`
• `` -> `...15`

performance <- read_excel("PilotPerformance.xlsx", sheet = 1)

New names:
• `` -> `...12`
• `` -> `...13`
• `` -> `...14`
• `` -> `...15`
• `` -> `...17`
• `` -> `...18`
• `` -> `...19`
• `` -> `...20`
• `` -> `...21`
• `` -> `...22`
• `` -> `...24`
• `` -> `...25`
• `` -> `...26`

cnvAmps <- read_csv("Data/ERPs/cnvAmps.csv")

Rows: 9 Columns: 5
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl (5): SubjectID, LoCog_LoMot, HiCog_LoMot, LoCog_HiMot, HiCog_HiMot

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

p3bAmps <- read_csv("Data/ERPs/p3bAmps.csv")

Rows: 9 Columns: 5
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl (5): SubjectID, LoCog_LoMot, HiCog_LoMot, LoCog_HiMot, HiCog_HiMot

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Data preparation

#Exclude Participants 9, 10, 18 due to technical issues (force plates, harness pain)
demographics <- demographics %>% filter(!ID %in% c(9, 10, 18))
demographics <- demographics %>% slice(1:9) # removing comments

nasa <- nasa %>% filter(!ID %in% c(9, 10, 18))
nasa <- nasa %>% select(1:12) # removing remarks

performance <- performance %>% filter(!ID %in% c(9, 10, 18))
performance <- performance %>% select(1:11) # removing comments

# checking data types
str(demographics)

tibble [9 × 43] (S3: tbl_df/tbl/data.frame)
 $ ID                : chr [1:9] "11" "12" "13" "14" ...
 $ Gender            : chr [1:9] "w" "m" "w" "m" ...
 $ Age               : chr [1:9] "25" "24" "22" "22" ...
 $ Height            : chr [1:9] "162" "173" "170" "170" ...
 $ Foot_width        : chr [1:9] "11" "12" "11" "12" ...
 $ Foot_length       : num [1:9] 27 32 28 29 29 28 29 34 28
 $ excluded          : chr [1:9] "0" "0" "0" "0" ...
 $ handedness        : chr [1:9] "right" "right" "right" "right" ...
 $ Vision Problems   : chr [1:9] "0" "0" "0" "0" ...
 $ MSSQ_Car          : num [1:9] 3 1 1 1 3 1 1 3 2
 $ MSSQ_Bus          : num [1:9] 2 1 1 0 4 2 2 3 1
 $ MSSQ_Train        : num [1:9] 1 2 1 0 2 1 1 3 1
 $ MSSQ_Plane        : num [1:9] 1 2 1 0 4 0 1 4 1
 $ MSSQ_Boat         : num [1:9] 0 1 1 0 1 3 1 2 0
 $ MSSQ_Ferry        : num [1:9] 0 3 1 0 2 3 1 1 0
 $ MSSQ_Swing        : num [1:9] 1 1 1 0 1 1 2 1 1
 $ MSSQ_Playground   : num [1:9] 1 1 1 2 2 1 1 1 3
 $ MSSQ_Rollercoaster: num [1:9] 1 1 1 1 0 2 0 1 1
 $ SS_Car            : num [1:9] 2 1 1 0 2 1 1 2 2
 $ SS_Bus            : num [1:9] 2 1 1 0 3 2 1 1 2
 $ SS_Train          : num [1:9] 1 2 1 0 1 1 1 2 2
 $ SS_Plane          : num [1:9] 1 3 1 0 1 1 1 1 1
 $ SS_Boat           : num [1:9] 0 3 2 0 1 2 1 2 0
 $ SS_Ferry          : num [1:9] 0 3 2 0 2 2 1 1 0
 $ SS_Swing          : num [1:9] 1 0 1 1 1 1 1 0 1
 $ SS_Playground     : num [1:9] 1 0 1 2 1 1 1 0 2
 $ SS_Rollercoaster  : num [1:9] 1 1 2 1 2 2 0 2 1
 $ SSQ_Overall       : num [1:9] 0 0 0 0 0 0 0 0 0
 $ SSQ_Tired         : num [1:9] 1 1 1 0 0 1 0 1 1
 $ SSQ_Headache      : num [1:9] 0 0 1 0 0 0 0 0 0
 $ SSQ_Eyetstrain    : num [1:9] 0 1 0 0 0 0 0 0 0
 $ SSQ_Sharpvision   : num [1:9] 0 2 0 0 0 0 0 0 0
 $ SSQ_Salvating     : num [1:9] 0 0 0 0 0 0 0 0 0
 $ SSQ_Sweating      : num [1:9] 0 0 0 1 1 0 0 0 0
 $ SSQ_Nauseous      : num [1:9] 0 0 0 0 0 0 0 0 0
 $ SSQ_Concentration : num [1:9] 0 1 0 0 0 0 0 0 0
 $ SSQ_Pressure      : num [1:9] 0 1 0 0 0 0 0 0 0
 $ SSQ_BlurredVision : num [1:9] 0 1 0 0 0 0 0 0 0
 $ SSQ_dizzy_open    : num [1:9] 0 0 0 0 0 0 0 0 0
 $ SSQ_dizzy_closed  : num [1:9] 0 0 0 0 0 0 0 0 0
 $ SSQ_Vertigo       : num [1:9] 0 0 0 0 0 0 0 0 0
 $ SSQ_Stomach       : num [1:9] 0 0 0 0 0 0 0 0 0
 $ SSQ_Burping       : num [1:9] 0 0 0 0 0 0 0 0 0

demographics <- demographics %>% 
  mutate(across(c(ID, Gender, `Vision Problems`), as.factor))
demographics <- demographics %>% 
  mutate(across(c(Age, Height, Foot_width, Foot_length), as.numeric))
str(nasa)

tibble [72 × 12] (S3: tbl_df/tbl/data.frame)
 $ ID         : num [1:72] 11 11 11 11 11 11 11 11 12 12 ...
 $ GA         : num [1:72] 14 16 11 15 12 10 8 10 17 17 ...
 $ KA         : num [1:72] 1 9 3 11 7 5 6 5 3 9 ...
 $ ZA         : num [1:72] 10 12 6 15 5 9 8 9 6 3 ...
 $ L          : num [1:72] 13 16 10 14 11 11 15 10 11 9 ...
 $ A          : num [1:72] 12 16 15 17 11 11 10 12 11 7 ...
 $ F          : num [1:72] 10 10 12 17 9 13 9 7 10 12 ...
 $ AG         : num [1:72] 1 5 2 3 6 2 2 3 2 3 ...
 $ rew_cog    : num [1:72] 1 1 2 2 1 1 2 2 1 1 ...
 $ rew_mot    : num [1:72] 1 2 1 2 1 2 1 2 1 2 ...
 $ RoundNumber: num [1:72] 1 1 1 1 2 2 2 2 1 1 ...
 $ BlockNumber: num [1:72] 3 4 1 2 6 5 8 7 2 3 ...

nasa <- nasa %>% 
  mutate(across(c(ID, rew_cog, rew_mot, RoundNumber, BlockNumber), as.factor))
str(performance)

tibble [81 × 11] (S3: tbl_df/tbl/data.frame)
 $ ID         : num [1:81] 11 11 11 11 11 11 11 11 11 12 ...
 $ RoundNumber: num [1:81] 0 1 1 1 1 2 2 2 2 0 ...
 $ rew_cog    : num [1:81] 0 1 1 2 2 1 1 2 2 0 ...
 $ rew_mot    : num [1:81] 0 1 2 1 2 1 2 1 2 0 ...
 $ cog_correct: num [1:81] 75.7 74 62.7 71.6 72.2 ...
 $ mot_correct: num [1:81] 91.8 89.9 90.3 91.4 90.7 ...
 $ L_I/A      : num [1:81] -5 -7 -7 -6 -6 -8 -9 -10 -9 -8 ...
 $ L_H/V      : num [1:81] 11 13 11 12 11 12 10 8 7 15 ...
 $ R_I/A      : num [1:81] 8 12 12 9 11 10 10 8 7 4 ...
 $ R_H/V      : num [1:81] 8 9 9 8 9 10 8 6 6 13 ...
 $ BlockNumber: num [1:81] 0 3 4 1 2 6 5 8 7 0 ...

performance <- performance %>% 
  mutate(across(c(ID, rew_cog, rew_mot, RoundNumber, BlockNumber), as.factor))
# removing baselines
performance <- performance %>% filter(RoundNumber != 0)

# EEG CNV + P3b
# transform into long format and duplicate rows to merge with data
colnames(cnvAmps)

[1] "SubjectID"   "LoCog_LoMot" "HiCog_LoMot" "LoCog_HiMot" "HiCog_HiMot"

cnv_long <- cnvAmps %>%  
  pivot_longer(
    cols = c(LoCog_LoMot, HiCog_LoMot, LoCog_HiMot, HiCog_HiMot), names_to = "reward_conditions", values_to = "cnv")
p3b_long <- p3bAmps %>%  
  pivot_longer(
    cols = c(LoCog_LoMot, HiCog_LoMot, LoCog_HiMot, HiCog_HiMot), names_to = "reward_conditions", values_to = "p3b")
# merging EEG dataset
eeg <- left_join(cnv_long, p3b_long, by = c("SubjectID", "reward_conditions"))

eeg <- eeg %>% 
  rename(
    ID = SubjectID
  )

eeg <- eeg %>% 
  mutate(reward_conditions = dplyr::recode(
    reward_conditions,
    "LoCog_LoMot" = "1_1",
    "HiCog_LoMot" = "2_1", 
    "LoCog_HiMot" = "1_2",
    "HiCog_HiMot" = "2_2"
  ))

# duplicating every entry in eeg to match data for merging
eeg <- eeg %>% 
  mutate(RoundNumber = 1) %>% 
  bind_rows(
    eeg %>%  mutate(RoundNumber = 2)
  )

str(eeg)

tibble [72 × 5] (S3: tbl_df/tbl/data.frame)
 $ ID               : num [1:72] 11 11 11 11 12 12 12 12 13 13 ...
 $ reward_conditions: chr [1:72] "1_1" "2_1" "1_2" "2_2" ...
 $ cnv              : num [1:72] -0.2467 -0.2522 -0.0856 -0.0707 -1.8224 ...
 $ p3b              : num [1:72] 0.0782 -0.2512 -0.4646 -0.3077 0.7029 ...
 $ RoundNumber      : num [1:72] 1 1 1 1 1 1 1 1 1 1 ...

eeg <- eeg %>% 
  mutate(across(c(ID, reward_conditions, RoundNumber), as.factor))

# joining NASA and performance scores by round/conditions
data <- left_join(performance, nasa, by = c("ID", "BlockNumber", "rew_cog", "rew_mot", "RoundNumber"))

data <- data %>% 
  mutate(reward_conditions = interaction(rew_cog, rew_mot, sep = "_"))

# joining eeg data and dataset (NASA + Performance)
data <- left_join(data, eeg, by = c("ID", "reward_conditions", "RoundNumber"))

Evaluating Questionnaires MSSQ, SSQ, NASA TLX

Motion Sickness Questionnaire MSSQ

# count number of methods of transport not used for MSSQ and SS questions (occurrences of 0)

demographics <- demographics %>% 
  mutate(
    MSSQ_not = rowSums(across(starts_with("MSSQ_"), ~ .x == 0)),
    SS_not = rowSums(across(starts_with("SS_"), ~ .x == 0))
  )

# transform MSSQ/SS to match scoring template (total remaining scored ones and subtract 9 to match scoring template)
demographics <- demographics %>%
  mutate(
    MSSQ_scored = rowSums(
      across(
        .cols = starts_with("MSSQ_") & !any_of("MSSQ_not")
      )
    ),
    SS_scored = rowSums(
      across(
        .cols = starts_with("SS_") & !any_of("SS_not")
      ) 
    )
  )

demographics <- demographics %>% 
  mutate(MSSQ_scored = MSSQ_scored - (9 - MSSQ_not),
         SS_scored = SS_scored - (9 - SS_not))

# formula
#MSSQ_child = total sickness score child * 9 / 9 * types not experienced
demographics <- demographics %>% 
  mutate(
    MSSQ_child = MSSQ_scored * 9 / 9 * MSSQ_not,
    MSSQ_adult = SS_scored * 9 / 9 * SS_not,
    MSSQ_total = MSSQ_child + MSSQ_adult
  )

# highest motion sickness score was a 14 (bottom 60 Percentile)
summary(demographics$MSSQ_total)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.000   6.000  10.000   8.556  12.000  14.000 

sd(demographics$MSSQ_total)

[1] 5.027701

Simulator Sickness Questionnaire SSQ

# gar nicht 0, stark = 3 
demographics <- demographics %>% 
  mutate(
    SSQ_N = rowSums(across(c(SSQ_Overall, SSQ_Salvating, SSQ_Sweating, SSQ_Nauseous, SSQ_Concentration, 
                             SSQ_Stomach, SSQ_Burping))) * 9.54,
    SSQ_O = rowSums(across(c(SSQ_Overall, SSQ_Tired, SSQ_Headache, SSQ_Eyetstrain, SSQ_Sharpvision, 
                             SSQ_Concentration, SSQ_BlurredVision))) * 7.58,
    SSQ_D = rowSums(across(c(SSQ_Sharpvision, SSQ_Nauseous, SSQ_Pressure, SSQ_BlurredVision, 
                             SSQ_dizzy_open, SSQ_dizzy_closed, SSQ_Vertigo))) * 13.92,
    SSQ_total = (SSQ_N + SSQ_O + SSQ_D) * 3.74
  )

summary(demographics$SSQ_total)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   0.00   28.35   28.35   72.83   35.68  414.02 

sd(demographics$SSQ_total)

[1] 128.7673

Subjective Workload NASA TLX

NASA TLX Descriptive Data

# NASA ratings per condition
data$mean_nasa <- rowMeans(data[, c("GA", "KA", "ZA", "L", "A", "F")])
    
mean_data <- data %>% 
  group_by(reward_conditions) %>% 
  summarise(
    mean_nasa_score = mean(mean_nasa, na.rm = TRUE),
    sd_nasa   = sd(mean_nasa, na.rm = TRUE),
    .groups = "drop"
  )
# Low/Low mean: 10.59259 sd: 3.131759
# High/Low mean: 11.75000 sd: 2.775506
# Low/High mean: 11.12963 sd: 3.351692
# High/High mean: 11.55556 sd: 3.585100

# Means and SDs of NASA Subscales
summary(data$GA) # mental demand

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   5.00   10.75   14.00   13.40   17.00   20.00 

sd(data$GA)

[1] 3.924139

summary(data$KA) # physical demand

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  1.000   5.750   9.500   9.431  13.000  20.000 

sd(data$KA)

[1] 5.026203

summary(data$ZA) # temporal demand

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   3.00    7.00   10.00   10.67   12.25   19.00 

sd(data$ZA)

[1] 4.275709

summary(data$L) # performance

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   5.00    9.00   12.00   12.04   15.25   20.00 

sd(data$L)

[1] 3.832542

summary(data$A) # effort

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    3.0    10.0    12.0    12.4    15.0    19.0 

sd(data$A)

[1] 3.71013

summary(data$F) # frustration

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  1.000   7.000   9.000   9.597  12.250  18.000 

sd(data$F)

[1] 4.297591

summary(data$AG) # gait complexity

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  1.000   3.000   8.500   8.653  13.000  19.000 

sd(data$AG)

[1] 5.828922

NASA TLX Linear Mixed Model + Interpretation

# Descriptive NASA TLX data
data <- data %>% select(ID, rew_cog, rew_mot, cog_correct, mot_correct, GA, KA, ZA, L, A, F, AG, reward_conditions, cnv, p3b)

str(data)

tibble [72 × 15] (S3: tbl_df/tbl/data.frame)
 $ ID               : Factor w/ 9 levels "11","12","13",..: 1 1 1 1 1 1 1 1 2 2 ...
 $ rew_cog          : Factor w/ 3 levels "0","1","2": 2 2 3 3 2 2 3 3 2 2 ...
 $ rew_mot          : Factor w/ 3 levels "0","1","2": 2 3 2 3 2 3 2 3 2 3 ...
 $ cog_correct      : num [1:72] 74 62.7 71.6 72.2 67.5 ...
 $ mot_correct      : num [1:72] 89.9 90.3 91.4 90.7 90 ...
 $ GA               : num [1:72] 14 16 11 15 12 10 8 10 17 17 ...
 $ KA               : num [1:72] 1 9 3 11 7 5 6 5 3 9 ...
 $ ZA               : num [1:72] 10 12 6 15 5 9 8 9 6 3 ...
 $ L                : num [1:72] 13 16 10 14 11 11 15 10 11 9 ...
 $ A                : num [1:72] 12 16 15 17 11 11 10 12 11 7 ...
 $ F                : num [1:72] 10 10 12 17 9 13 9 7 10 12 ...
 $ AG               : num [1:72] 1 5 2 3 6 2 2 3 2 3 ...
 $ reward_conditions: Factor w/ 9 levels "0_0","1_0","2_0",..: 5 8 6 9 5 8 6 9 5 8 ...
 $ cnv              : num [1:72] -0.2467 -0.0856 -0.2522 -0.0707 -0.2467 ...
 $ p3b              : num [1:72] 0.0782 -0.4646 -0.2512 -0.3077 0.0782 ...

data$rew_cog <- factor(data$rew_cog, levels = c(1, 2), labels = c("Low", "High"))
data$rew_mot <- factor(data$rew_mot, levels = c(1, 2), labels = c("Low", "High"))
data$reward_conditions <- interaction(data$rew_cog, data$rew_mot, sep = "/")

# separate model for each NASA TLX and Gait Complexity subscale

# mental demand, physical demand, temporal demand, performance, effort, frustration, gait complexity
scales <- c("GA", "KA", "ZA", "L", "A", "F", "AG")

# All models follow the same structure: 
# DV: score on that scale, IV: RewardLevelStock, RewardLevelWalking, Random factor participant ID

nasa_models <- map(scales, ~ lmer(as.formula(paste0(.x, " ~ rew_cog * rew_mot + (1|ID)")), 
                                  data = data)) |> set_names(scales)
anova_all <- imap_dfr(nasa_models, ~ {
  Anova(.x, type = 3) |>
    as.data.frame() |>
    rownames_to_column("Effect") |>
    mutate(Scale = .y)
})

# keep reward terms, format nicely
order_effects <- c("rew_cog","rew_mot","rew_cog:rew_mot")
fmt_p <- function(p) ifelse(p < .001, "< .001", sprintf("%.3f", p))
stars <- function(p) case_when(p < .001 ~ "***",
                               p < .01  ~ "**",
                               p < .05  ~ "*",
                               TRUE     ~ "")

anova_view <- anova_all %>%
  filter(Effect %in% order_effects) %>%
  mutate(Effect = factor(Effect, levels = order_effects),
         ChiSq  = round(Chisq, 2),
         p      = fmt_p(`Pr(>Chisq)`),
         Sig    = stars(`Pr(>Chisq)`)) %>%
  arrange(Scale, Effect) %>%
  select(Scale, Effect, df = Df, ChiSq, p, Sig)

# print a clean table
kable(anova_view, align = "l")

Scale	Effect	df	ChiSq	p	Sig
A	rew_cog	1	4.86	0.027	*
A	rew_mot	1	0.09	0.766
A	rew_cog:rew_mot	1	0.21	0.643
AG	rew_cog	1	3.36	0.067
AG	rew_mot	1	1.94	0.164
AG	rew_cog:rew_mot	1	0.60	0.437
F	rew_cog	1	1.43	0.231
F	rew_mot	1	0.94	0.332
F	rew_cog:rew_mot	1	0.27	0.600
GA	rew_cog	1	5.20	0.023	*
GA	rew_mot	1	0.10	0.751
GA	rew_cog:rew_mot	1	0.58	0.446
KA	rew_cog	1	0.00	1.000
KA	rew_mot	1	4.25	0.039	*
KA	rew_cog:rew_mot	1	1.03	0.311
L	rew_cog	1	1.80	0.179
L	rew_mot	1	0.01	0.941
L	rew_cog:rew_mot	1	2.03	0.154
ZA	rew_cog	1	1.27	0.261
ZA	rew_mot	1	0.27	0.604
ZA	rew_cog:rew_mot	1	0.73	0.392

nasa_summary_view <- imap_dfr(nasa_models, ~{
  sm <- summary(.x)
  tb <- as.data.frame(coef(sm))
  tb$Term <- rownames(tb)
  rownames(tb) <- NULL
  tb$Scale <- .y
  tb
}) %>%
  relocate(Scale, Term) %>%
  rename(
    Estimate = Estimate,
    SE       = `Std. Error`,
    df       = df,
    t        = `t value`,
    p_raw    = `Pr(>|t|)`
  ) %>%
  mutate(
    p  = if ("p_raw" %in% names(.)) fmt_p(p_raw) else NA_character_,
    Sig = if ("p_raw" %in% names(.)) stars(p_raw) else ""
  ) %>%
  select(Scale, Term, Estimate, SE, df, t, p, Sig)

kable(nasa_summary_view, align = "l")

Scale	Term	Estimate	SE	df	t	p	Sig
GA	(Intercept)	12.5000000	1.1765660	12.649665	10.6241386	< .001	***
GA	rew_cogHigh	2.0000000	0.8770919	60.000000	2.2802629	0.026	*
GA	rew_motHigh	0.2777778	0.8770919	60.000000	0.3167032	0.753
GA	rew_cogHigh:rew_motHigh	-0.9444444	1.2403952	60.000000	-0.7614061	0.449
KA	(Intercept)	8.7777778	1.5816390	10.802449	5.5497987	< .001	***
KA	rew_cogHigh	0.0000000	0.9697588	60.000000	0.0000000	1.000
KA	rew_motHigh	2.0000000	0.9697588	60.000000	2.0623685	0.044	*
KA	rew_cogHigh:rew_motHigh	-1.3888889	1.3714460	60.000000	-1.0127186	0.315
ZA	(Intercept)	10.2777778	1.4065201	9.403966	7.3072386	< .001	***
ZA	rew_cogHigh	0.7222222	0.6419011	60.000000	1.1251300	0.265
ZA	rew_motHigh	-0.3333333	0.6419011	60.000000	-0.5192908	0.605
ZA	rew_cogHigh:rew_motHigh	0.7777778	0.9077853	60.000000	0.8567861	0.395
L	(Intercept)	11.8888889	1.2100762	10.825006	9.8249094	< .001	***
L	rew_cogHigh	1.0000000	0.7443546	60.000000	1.3434457	0.184
L	rew_motHigh	0.0555556	0.7443546	60.000000	0.0746359	0.941
L	rew_cogHigh:rew_motHigh	-1.5000000	1.0526764	60.000000	-1.4249393	0.159
A	(Intercept)	11.3888889	1.0533237	15.670583	10.8123353	< .001	***
A	rew_cogHigh	2.0555556	0.9319986	60.000000	2.2055350	0.031	*
A	rew_motHigh	0.2777778	0.9319986	60.000000	0.2980453	0.767
A	rew_cogHigh:rew_motHigh	-0.6111111	1.3180451	60.000000	-0.4636496	0.645
F	(Intercept)	8.7222222	1.3069216	12.645387	6.6738681	< .001	***
F	rew_cogHigh	1.1666667	0.9739403	60.000000	1.1978832	0.236
F	rew_motHigh	0.9444444	0.9739403	60.000000	0.9697149	0.336
F	rew_cogHigh:rew_motHigh	-0.7222222	1.3773596	60.000000	-0.5243527	0.602
AG	(Intercept)	8.6111111	1.9396419	8.996932	4.4395365	0.002	**
AG	rew_cogHigh	-1.3888889	0.7579483	60.000000	-1.8324323	0.072
AG	rew_motHigh	1.0555556	0.7579483	60.000000	1.3926485	0.169
AG	rew_cogHigh:rew_motHigh	0.8333333	1.0719007	60.000000	0.7774352	0.440

# get EMMs per model and combine
emmeans_nasa_models <- imap_dfr(nasa_models, ~ as.data.frame(emmeans(.x, ~ rew_cog * rew_mot)) %>%
                      mutate(Scale = .y))

# plot all effects together
nasa_ploting <- emmeans_nasa_models %>%
  mutate(cond = interaction(rew_cog, rew_mot, sep = "/"))

nasa_ploting$Scale <- factor(nasa_ploting$Scale,
                                   levels = c("GA", "KA", "ZA", "L", "A", "F", "AG"))
str(nasa_ploting)

'data.frame':   28 obs. of  9 variables:
 $ rew_cog : Factor w/ 2 levels "Low","High": 1 2 1 2 1 2 1 2 1 2 ...
 $ rew_mot : Factor w/ 2 levels "Low","High": 1 1 2 2 1 1 2 2 1 1 ...
 $ emmean  : num  12.5 14.5 12.78 13.83 8.78 ...
 $ SE      : num  1.18 1.18 1.18 1.18 1.58 ...
 $ df      : num  12.6 12.6 12.6 12.6 10.8 ...
 $ lower.CL: num  9.95 11.95 10.23 11.28 5.29 ...
 $ upper.CL: num  15 17 15.3 16.4 12.3 ...
 $ Scale   : Factor w/ 7 levels "GA","KA","ZA",..: 1 1 1 1 2 2 2 2 3 3 ...
 $ cond    : Factor w/ 4 levels "Low/Low","High/Low",..: 1 2 3 4 1 2 3 4 1 2 ...

nasa_effects_plot <- ggplot(nasa_ploting, aes(x = Scale, y = emmean,
               color = cond, shape = cond, group = cond)) +
  geom_point(position = position_dodge(width = 0.35), size = 2) +
  geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL),
                position = position_dodge(width = 0.35), width = 0.2) +
  labs(
    y = "Average Subscale scores (1–20)",
    x = "NASA TLX and Gait Complexity Subscales",
    color = "Reward condition\n (Stock/Walking)",
    shape = "Reward condition\n (Stock/Walking)"
  ) +
  scale_shape_manual(values = c("Low/Low" = 16, "Low/High" = 16,
                                "High/Low" = 17, "High/High" = 17)) +
  scale_color_manual(values = c("Low/Low" = "#1f77b4",   # motor Low
                                "High/Low" = "#1f77b4",  # motor Low
                                "Low/High" = "#d62728",  # motor High
                                "High/High" = "#d62728")) + # motor High
  scale_x_discrete(labels = c("GA" = "Mental Demand", "KA" = "Physical Demand", "ZA" = "Temporal Demand", "L" = "Performance", "A" = "Effort", "F" = "Frustration", "AG" = "Gait Complexity"))

ggsave("../../ReStoWa/MasterThesis/Graphs/nasa_tlx_conditions.jpg", nasa_effects_plot, width = 11, height = 4, dpi = 300)

Performance Data Descriptives

# editing to make variables easier to understand
data <- data %>% 
  rename(
    WalkingAccuracy = mot_correct,
    StockAccuracy = cog_correct
  )

# descriptive performance data

mean_perf <- data %>%
  group_by(reward_conditions) %>%
  summarise(
    mean_stock   = mean(StockAccuracy, na.rm = TRUE),
    mean_walking = mean(WalkingAccuracy, na.rm = TRUE),
    sd_walk = sd(WalkingAccuracy),
    sd_stock = sd(StockAccuracy),
    .groups = "drop"
  )

# Order: LowLow LowHigh HighLow HighHigh
# Stock 68.21278, 70.28222, 66.66889, 71.56500 (in order: HH, LH, LL, HL)
# Walk: 87.00389, 87.98333, 88.38444, 71.56500 (high to low in order: HL, LH, LL, HH)

Stock Monitoring Task Performance Linear Mixed Model

# DV: StockAccuracy, IV: StockRewardLevel, WalkingRewardLevel, Random factor: participant ID

perf_stock_model <- lmer(StockAccuracy ~ rew_cog * rew_mot + (1|ID), data = data, REML = FALSE)

Anova(perf_stock_model) # Stock task significant

Analysis of Deviance Table (Type II Wald chisquare tests)

Response: StockAccuracy
                 Chisq Df Pr(>Chisq)   
rew_cog         6.8596  1   0.008817 **
rew_mot         0.0096  1   0.921790   
rew_cog:rew_mot 1.1296  1   0.287853   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(perf_stock_model)

Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
  method [lmerModLmerTest]
Formula: StockAccuracy ~ rew_cog * rew_mot + (1 | ID)
   Data: data

      AIC       BIC    logLik -2*log(L)  df.resid 
    491.8     505.5    -239.9     479.8        66 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.6357 -0.4791  0.0311  0.6221  1.5559 

Random effects:
 Groups   Name        Variance Std.Dev.
 ID       (Intercept) 70.35    8.388   
 Residual             31.83    5.642   
Number of obs: 72, groups:  ID, 9

Fixed effects:
                        Estimate Std. Error     df t value Pr(>|t|)    
(Intercept)               68.213      3.096 12.078  22.032 4.03e-11 ***
rew_cogHigh                2.069      1.881 63.000   1.100    0.275    
rew_motHigh               -1.544      1.881 63.000  -0.821    0.415    
rew_cogHigh:rew_motHigh    2.827      2.660 63.000   1.063    0.292    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) rw_cgH rw_mtH
rew_cogHigh -0.304              
rew_motHigh -0.304  0.500       
rw_cgHgh:_H  0.215 -0.707 -0.707

perf_stock_emmeans <- emmeans(perf_stock_model, pairwise ~ rew_cog * rew_mot)
print(perf_stock_emmeans) # only normal - short significant

$emmeans
 rew_cog rew_mot emmean   SE   df lower.CL upper.CL
 Low     Low       68.2 3.27 13.6     61.2     75.2
 High    Low       70.3 3.27 13.6     63.3     77.3
 Low     High      66.7 3.27 13.6     59.6     73.7
 High    High      71.6 3.27 13.6     64.5     78.6

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

$contrasts
 contrast             estimate   SE   df t.ratio p.value
 Low Low - High Low      -2.07 1.93 66.2  -1.074  0.7065
 Low Low - Low High       1.54 1.93 66.2   0.801  0.8536
 Low Low - High High     -3.35 1.93 66.2  -1.740  0.3118
 High Low - Low High      3.61 1.93 66.2   1.875  0.2486
 High Low - High High    -1.28 1.93 66.2  -0.666  0.9096
 Low High - High High    -4.90 1.93 66.2  -2.541  0.0628

Degrees-of-freedom method: kenward-roger 
P value adjustment: tukey method for comparing a family of 4 estimates 

perf_stock_effects <- allEffects(perf_stock_model)
perf_stock_effects_model <- as.data.frame(perf_stock_effects[[1]])

perf_stock_plot <- 
  ggplot(perf_stock_effects_model,
       aes(x = rew_mot, y = fit,
           color = rew_cog, group = rew_cog)) +
  geom_point(position = position_dodge(width = 0.35), size = 3) +
  geom_errorbar(aes(ymin = fit - se, ymax = fit + se),
                position = position_dodge(width = 0.35), width = 0.3, size = 1) +
  ylab("Mean Stock Monitoring Task Accuracy in (%)") +
  xlab("Walking Task Reward Level") +
  labs(color = "Stock Monitoring Task\nReward Level") + 
  scale_color_discrete(labels = c("Low", "High")
  )

Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.

plot(perf_stock_plot)

ggsave("../../ReStoWa/MasterThesis/Graphs/perf_stock_plot.jpg", perf_stock_plot, width = 11, height = 4, dpi = 300)

Walking Task linear mixed model

# DV: WalkingAccuracy, IV: StockRewardLevel, WalkingRewardLevel, Random factor: participant ID

perf_walk_model <- lmer(WalkingAccuracy ~ rew_cog * rew_mot + (1|ID), data = data, REML = FALSE)

Anova(perf_walk_model) # nothing significant

Analysis of Deviance Table (Type II Wald chisquare tests)

Response: WalkingAccuracy
                 Chisq Df Pr(>Chisq)
rew_cog         0.5558  1     0.4560
rew_mot         1.1124  1     0.2916
rew_cog:rew_mot 0.0001  1     0.9925

summary(perf_walk_model)

Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
  method [lmerModLmerTest]
Formula: WalkingAccuracy ~ rew_cog * rew_mot + (1 | ID)
   Data: data

      AIC       BIC    logLik -2*log(L)  df.resid 
    483.3     497.0    -235.6     471.3        66 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.6203 -0.1356  0.1253  0.3137  2.2415 

Random effects:
 Groups   Name        Variance Std.Dev.
 ID       (Intercept) 36.9     6.075   
 Residual             30.3     5.504   
Number of obs: 72, groups:  ID, 9

Fixed effects:
                        Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)             87.00389    2.40482 14.56572  36.179 1.13e-15 ***
rew_cogHigh              0.97944    1.83475 63.00000   0.534    0.595    
rew_motHigh              1.38056    1.83475 63.00000   0.752    0.455    
rew_cogHigh:rew_motHigh -0.02444    2.59473 63.00000  -0.009    0.993    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) rw_cgH rw_mtH
rew_cogHigh -0.381              
rew_motHigh -0.381  0.500       
rw_cgHgh:_H  0.270 -0.707 -0.707

perf_walk_emmeans <- emmeans(perf_walk_model, pairwise ~ rew_cog * rew_mot)
print(perf_walk_emmeans) # no significant results here

$emmeans
 rew_cog rew_mot emmean   SE   df lower.CL upper.CL
 Low     Low       87.0 2.53 16.4     81.6     92.4
 High    Low       88.0 2.53 16.4     82.6     93.3
 Low     High      88.4 2.53 16.4     83.0     93.7
 High    High      89.3 2.53 16.4     84.0     94.7

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

$contrasts
 contrast             estimate   SE   df t.ratio p.value
 Low Low - High Low     -0.979 1.88 66.2  -0.521  0.9538
 Low Low - Low High     -1.381 1.88 66.2  -0.734  0.8829
 Low Low - High High    -2.336 1.88 66.2  -1.242  0.6025
 High Low - Low High    -0.401 1.88 66.2  -0.213  0.9965
 High Low - High High   -1.356 1.88 66.2  -0.721  0.8883
 Low High - High High   -0.955 1.88 66.2  -0.508  0.9569

Degrees-of-freedom method: kenward-roger 
P value adjustment: tukey method for comparing a family of 4 estimates 

perf_walk_effects <- allEffects(perf_walk_model)
perf_walk_effects_model <- as.data.frame(perf_walk_effects[[1]])

perf_walk_plot <- 
  ggplot(perf_walk_effects_model,
       aes(x = rew_mot, y = fit,
           color = rew_cog, group = rew_cog)) +
  geom_point(position = position_dodge(width = 0.35), size = 3) +
  geom_errorbar(aes(ymin = fit - se, ymax = fit + se),
                position = position_dodge(width = 0.35), width = 0.3, size = 1) +
  ylab("Mean Walking Task Accuracy in (%)") +
  xlab("Walking Task Reward Level") +
  labs(color = "Stock Monitoring Task\nReward Level") + 
  scale_color_discrete(labels = c("Low", "High")
  )

plot(perf_walk_plot)

ggsave("../../ReStoWa/MasterThesis/Graphs/perf_walk_plot.jpg", perf_walk_plot, width = 11, height = 4, dpi = 300)

EEG Data descriptives

# mean amplitude per condition 

mean_eeg_components <- data %>%
  group_by(reward_conditions) %>%
  summarise(
    mean_cnv = mean(cnv, na.rm = TRUE),
    sd_cnv   = sd(cnv, na.rm = TRUE),
    mean_p3b = mean(p3b, na.rm = TRUE),
    sd_p3b   = sd(p3b, na.rm = TRUE),
    .groups = "drop"
  )

# In order: LowLow HighLow HighLow HighHigh
# CNV: -0.8130322 -1.0749466 -0.7426372 -1.1539051 (Most negative)
  # More negative = stronger preparation (best to worst)
  # HH, HL, LL, LH
# p3b: 0.5699609 0.2949174 0.2128196 0.1913472
  # more positive = stronger attention/processing (best to worst)
  # LL, HL, LH, HH

P3b linear mixed model analysis

p3b_model <- lmer(p3b ~ rew_cog * rew_mot + (1|ID), data = data, REML = FALSE)

Anova(p3b_model) # nothing significant

Analysis of Deviance Table (Type II Wald chisquare tests)

Response: p3b
                 Chisq Df Pr(>Chisq)
rew_cog         0.6862  1     0.4075
rew_mot         1.6565  1     0.1981
rew_cog:rew_mot 0.5018  1     0.4787

summary(p3b_model)

Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
  method [lmerModLmerTest]
Formula: p3b ~ rew_cog * rew_mot + (1 | ID)
   Data: data

      AIC       BIC    logLik -2*log(L)  df.resid 
    198.8     212.4     -93.4     186.8        66 

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.88977 -0.38305  0.01252  0.46111  2.85783 

Random effects:
 Groups   Name        Variance Std.Dev.
 ID       (Intercept) 0.7653   0.8748  
 Residual             0.5766   0.7593  
Number of obs: 72, groups:  ID, 9

Fixed effects:
                        Estimate Std. Error      df t value Pr(>|t|)
(Intercept)               0.5700     0.3422 14.1134   1.666    0.118
rew_cogHigh              -0.2750     0.2531 63.0000  -1.087    0.281
rew_motHigh              -0.3571     0.2531 63.0000  -1.411    0.163
rew_cogHigh:rew_motHigh   0.2536     0.3580 63.0000   0.708    0.481

Correlation of Fixed Effects:
            (Intr) rw_cgH rw_mtH
rew_cogHigh -0.370              
rew_motHigh -0.370  0.500       
rw_cgHgh:_H  0.262 -0.707 -0.707

p3b_emmeans <- emmeans(p3b_model, pairwise ~ rew_cog * rew_mot)
print(p3b_emmeans) # nothing significant 

$emmeans
 rew_cog rew_mot emmean   SE   df lower.CL upper.CL
 Low     Low      0.570 0.36 15.9   -0.195    1.335
 High    Low      0.295 0.36 15.9   -0.470    1.059
 Low     High     0.213 0.36 15.9   -0.552    0.977
 High    High     0.191 0.36 15.9   -0.573    0.956

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

$contrasts
 contrast             estimate    SE   df t.ratio p.value
 Low Low - High Low     0.2750 0.259 66.2   1.060  0.7146
 Low Low - Low High     0.3571 0.259 66.2   1.377  0.5181
 Low Low - High High    0.3786 0.259 66.2   1.460  0.4674
 High Low - Low High    0.0821 0.259 66.2   0.317  0.9889
 High Low - High High   0.1036 0.259 66.2   0.399  0.9783
 Low High - High High   0.0215 0.259 66.2   0.083  0.9998

Degrees-of-freedom method: kenward-roger 
P value adjustment: tukey method for comparing a family of 4 estimates 

p3b_effects <- allEffects(p3b_model)
p3b_effects_model <- as.data.frame(p3b_effects[[1]])

p3b_plot <- 
  ggplot(p3b_effects_model,
       aes(x = rew_mot, y = fit,
           color = rew_cog, group = rew_cog)) +
  geom_point(position = position_dodge(width = 0.35), size = 3) +
  geom_errorbar(aes(ymin = fit - se, ymax = fit + se),
                position = position_dodge(width = 0.35), width = 0.3, size = 1) +
  ylab("Mean P3b Amplitude in (µV)") +
  xlab("Walking Task Reward Level") +
  labs(color = "Stock Monitoring Task\nReward Level") + 
  scale_color_discrete(labels = c("Low", "High")
  )

plot(p3b_plot)

ggsave("../../ReStoWa/MasterThesis/Graphs/p3b_plot.jpg", p3b_plot, width = 11, height = 4, dpi = 300)

CNV linear mixed model analysis

cnv_model <- lmer(cnv ~ rew_cog * rew_mot + (1|ID), data = data, REML = FALSE)

Anova(cnv_model) # StockReward, nothing else

Analysis of Deviance Table (Type II Wald chisquare tests)

Response: cnv
                 Chisq Df Pr(>Chisq)  
rew_cog         5.6522  1    0.01743 *
rew_mot         0.0009  1    0.97587  
rew_cog:rew_mot 0.2782  1    0.59787  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(cnv_model)

Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
  method [lmerModLmerTest]
Formula: cnv ~ rew_cog * rew_mot + (1 | ID)
   Data: data

      AIC       BIC    logLik -2*log(L)  df.resid 
    173.1     186.8     -80.5     161.1        66 

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.65392 -0.32826  0.04036  0.48670  2.38656 

Random effects:
 Groups   Name        Variance Std.Dev.
 ID       (Intercept) 1.2434   1.1151  
 Residual             0.3608   0.6007  
Number of obs: 72, groups:  ID, 9

Fixed effects:
                        Estimate Std. Error      df t value Pr(>|t|)  
(Intercept)              -0.8130     0.3977 10.9720  -2.044   0.0657 .
rew_cogHigh              -0.2619     0.2002 63.0000  -1.308   0.1956  
rew_motHigh               0.0704     0.2002 63.0000   0.352   0.7263  
rew_cogHigh:rew_motHigh  -0.1494     0.2832 63.0000  -0.527   0.5997  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) rw_cgH rw_mtH
rew_cogHigh -0.252              
rew_motHigh -0.252  0.500       
rw_cgHgh:_H  0.178 -0.707 -0.707

cnv_emmeans <- emmeans(cnv_model, pairwise ~ rew_cog * rew_mot)
print(cnv_emmeans) # no significance, closest LH - HH

$emmeans
 rew_cog rew_mot emmean    SE   df lower.CL upper.CL
 Low     Low     -0.813 0.421 12.3    -1.73    0.100
 High    Low     -1.075 0.421 12.3    -1.99   -0.161
 Low     High    -0.743 0.421 12.3    -1.66    0.171
 High    High    -1.154 0.421 12.3    -2.07   -0.240

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

$contrasts
 contrast             estimate    SE   df t.ratio p.value
 Low Low - High Low     0.2619 0.205 66.2   1.277  0.5809
 Low Low - Low High    -0.0704 0.205 66.2  -0.343  0.9860
 Low Low - High High    0.3409 0.205 66.2   1.661  0.3521
 High Low - Low High   -0.3323 0.205 66.2  -1.620  0.3748
 High Low - High High   0.0790 0.205 66.2   0.385  0.9805
 Low High - High High   0.4113 0.205 66.2   2.005  0.1967

Degrees-of-freedom method: kenward-roger 
P value adjustment: tukey method for comparing a family of 4 estimates 

cnv_effects <- allEffects(cnv_model)
cnv_effects_model <- as.data.frame(cnv_effects[[1]])

cnv_plot <- 
  ggplot(cnv_effects_model,
       aes(x = rew_mot, y = fit,
           color = rew_cog, group = rew_cog)) +
  geom_point(position = position_dodge(width = 0.35), size = 3) +
  geom_errorbar(aes(ymin = fit - se, ymax = fit + se),
                position = position_dodge(width = 0.35), width = 0.3, size = 1) +
  ylab("Mean CNV Amplitude in (µV)") +
  xlab("Walking Task Reward Level") +
  labs(color = "Stock Monitoring Task\nReward Level") + 
  scale_color_discrete(labels = c("Low", "High")
  )

plot(cnv_plot)

ggsave("../../ReStoWa/MasterThesis/Graphs/cnv_plot.jpg", cnv_plot, width = 11, height = 4, dpi = 300)

Visualizing CNV grand average means

# Load
cnv_means <- read_csv("Data/ERPs/cnvMeansTable.csv")

Rows: 4 Columns: 3002
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr    (2): SplitLabel, ConditionLabel
dbl (3000): MeanData_1, MeanData_2, MeanData_3, MeanData_4, MeanData_5, Mean...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

cnv_times <- read_csv("Data/ERPs/cnvMeansTimes.csv")

Rows: 0 Columns: 1500
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1500): -1500, -1498, -1496, -1494, -1492, -1490, -1488, -1486, -1484, -...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

# Extract timepoints from colnames of the times file
timepoints <- as.numeric(colnames(cnv_times))

# Reshape mean amplitudes
cnv_long <- cnv_means %>%
  pivot_longer(cols = starts_with("MeanData_"),
               names_to = "time_index", values_to = "mean") %>%
  mutate(time = timepoints[as.numeric(str_remove(time_index, "MeanData_"))]) %>%
  select(SplitLabel, ConditionLabel, time, mean)

# Reshape SEM
cnv_sem <- cnv_means %>%
  pivot_longer(cols = starts_with("SemData_"),
               names_to = "time_index", values_to = "sem") %>%
  mutate(time = timepoints[as.numeric(str_remove(time_index, "SemData_"))]) %>%
  select(SplitLabel, ConditionLabel, time, sem)

# Merge
cnv_long <- left_join(cnv_long, cnv_sem,
                      by = c("SplitLabel", "ConditionLabel", "time"))

# Plot CNV

cnv_long <- cnv_long %>%
  mutate(
    SplitLabel = factor(SplitLabel, levels = c("loCog", "hiCog"), 
                        labels = c("Low Reward Stock Monitoring Task", "High Reward Stock Monitoring Task")),
    ConditionLabel = factor(ConditionLabel, levels = c("loMot", "hiMot"),
                            labels = c("Low Reward", "High Reward"))
  )

# if I understand it correctly we have the CNV peak at -100 (windows of 100ms: -200 to 0ms)

cnv_means_plot <- ggplot(cnv_long, aes(x = time, y = mean, color = ConditionLabel)) +
  geom_rect(
    data = data.frame(xmin = -200, xmax = 0, ymin = -Inf, ymax = Inf),
    aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax),
    inherit.aes = FALSE, fill = "red", alpha = 0.1
  ) +
  geom_line(size = 1) +
  geom_ribbon(aes(ymin = mean - sem, ymax = mean + sem, fill = ConditionLabel),
              alpha = 0.2, color = NA) +
  facet_wrap(~ SplitLabel, ncol = 1) +
  labs(x = "Time in (ms)", 
       y = "Grand Mean CNV Amplitudes in (µV)", 
       color = "Walking Task Reward Levels", 
       fill  = "Walking Task Reward Levels") +
  theme_minimal()

ggsave("Graphs/cnv_means_plot.jpg", cnv_means_plot, width = 11, height = 4, dpi = 300)

Visualizing P3b grand average means

p3_means <- read_csv("Data/ERPs/p3MeansTable.csv")

Rows: 4 Columns: 3002
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr    (2): SplitLabel, ConditionLabel
dbl (3000): MeanData_1, MeanData_2, MeanData_3, MeanData_4, MeanData_5, Mean...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

p3_times <- read_csv("Data/ERPs/p3MeansTimes.csv")

Rows: 0 Columns: 1500
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1500): -1500, -1498, -1496, -1494, -1492, -1490, -1488, -1486, -1484, -...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

timepoints <- as.numeric(colnames(p3_times))

# Long format for plotting
p3_long <- p3_means %>%
  pivot_longer(cols = starts_with("MeanData_"),
               names_to = "time_index", values_to = "mean") %>%
  mutate(time = timepoints[as.numeric(str_remove(time_index, "MeanData_"))]) %>%
  select(SplitLabel, ConditionLabel, time, mean)

# Reshape Standard erorr of the mean
p3_sem <- p3_means %>%
  pivot_longer(cols = starts_with("SemData_"),
               names_to = "time_index", values_to = "sem") %>%
  mutate(time = timepoints[as.numeric(str_remove(time_index, "SemData_"))]) %>%
  select(SplitLabel, ConditionLabel, time, sem)

p3_long <- left_join(p3_long, p3_sem,
                     by = c("SplitLabel", "ConditionLabel", "time"))

# Change factors for order within graphs
p3_long <- p3_long %>%
  mutate(
    SplitLabel = factor(SplitLabel, levels = c("loCog", "hiCog"), 
                        labels = c("Low Reward Stock Monitoring Task", "High Reward Stock Monitoring Task")),
    ConditionLabel = factor(ConditionLabel, levels = c("loMot", "hiMot"),
                            labels = c("Low Reward", "High Reward"))
  )

# P3b peak at 844 (windows of 100ms 744-944)

p3_means_plot <- ggplot(p3_long, aes(x = time, y = mean, color = ConditionLabel)) +
  geom_rect(
    data = data.frame(xmin = 744, xmax = 944, ymin = -Inf, ymax = Inf), # typical P3b window
    aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax),
    inherit.aes = FALSE, fill = "red", alpha = 0.1
  ) +
  geom_line(size = 1) +
  geom_ribbon(aes(ymin = mean - sem, ymax = mean + sem, fill = ConditionLabel),
              alpha = 0.2, color = NA) +
  facet_wrap(~ SplitLabel, ncol = 1) +
  labs(x = "Time in (ms)", 
       y = "Grand Mean P3b Amplitudes in (µV)", 
       color = "Walking Task Reward Levels", 
       fill  = "Walking Task Reward Levels") +
  theme_minimal()

ggsave("Graphs/p3_means_plot.jpg", p3_means_plot, width = 11, height = 4, dpi = 300)