Thesis Analysis

Loading Packages

library(tidyverse)

## ── 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.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.0.4     
## ── 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(readr)
library(ggplot2)
library(tidyr)
library(dplyr)
library(lme4)

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

library(emmeans)

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

library(rstanarm)

## Loading required package: Rcpp
## This is rstanarm version 2.32.1
## - See https://mc-stan.org/rstanarm/articles/priors for changes to default priors!
## - Default priors may change, so it's safest to specify priors, even if equivalent to the defaults.
## - For execution on a local, multicore CPU with excess RAM we recommend calling
##   options(mc.cores = parallel::detectCores())

library(ggplot2)
library(effects)

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

library(car)

## 
## Attaching package: 'car'
## 
## The following object is masked from 'package:rstanarm':
## 
##     logit
## 
## The following object is masked from 'package:dplyr':
## 
##     recode
## 
## The following object is masked from 'package:purrr':
## 
##     some

Importing Data

Experiment <-read.csv("all_excluded.csv", sep = ";" )
Questionnaire <-read.csv("Questionaire_Data3.csv", sep = ";")

Adjusting the Questionnaire Dataset

# Deleting Columns (Run codes seperately)
Questionnaire <- Questionnaire[, -c(1:4)]
Questionnaire <- Questionnaire[, -c(2:13)]

# Renaming Columns
colnames(Questionnaire)[colnames(Questionnaire) == "Handedness"] <- "Footedness"
colnames(Questionnaire)[colnames(Questionnaire) == "Q24"] <- "Gender"

# Removing Participants (To Merge with other Dataset)
Questionnaire <- Questionnaire %>% 
  filter(!ID %in% c(6, 9, 21))

# Adjusting weight of ID 8
Questionnaire$Weight[Questionnaire$ID == 8] <- 62

## Creating New Variables
# Pleasure_B and Arousal_B
Questionnaire$Pleasure_B <- as.numeric(sapply(strsplit(as.character(Questionnaire$C_B), "/"), `[`, 1))
Questionnaire$Arousal_B <- as.numeric(sapply(strsplit(as.character(Questionnaire$C_B), "/"), `[`, 2))


# Pleasure_1 and Arousal_1
Questionnaire$Pleasure_1 <- sapply(strsplit(as.character(Questionnaire$C_1), ","), function(x) {
  parts <- strsplit(x[1], "/")[[1]]
  return(as.numeric(parts[1])) 
})

Questionnaire$Arousal_1 <- sapply(strsplit(as.character(Questionnaire$C_1), ","), function(x) {
  parts <- strsplit(x[length(x)], "/")[[1]]
  return(as.numeric(parts[2]))  
})


# Pleasure_2 and Arousal_2
Questionnaire$Pleasure_2 <- as.numeric(sapply(strsplit(as.character(Questionnaire$C_2), "/"), `[`, 1))
Questionnaire$Arousal_2 <- as.numeric(sapply(strsplit(as.character(Questionnaire$C_2), "/"), `[`, 2))

# Pleasure_3 and Arousal_3
Questionnaire$Pleasure_3 <- as.numeric(sapply(strsplit(as.character(Questionnaire$C_3), "/"), `[`, 1))
Questionnaire$Arousal_3 <- as.numeric(sapply(strsplit(as.character(Questionnaire$C_3), "/"), `[`, 2))


# Pleasure_4 and Arousal_4
Questionnaire$Pleasure_4 <- sapply(strsplit(as.character(Questionnaire$C_4), ","), function(x) {
  parts <- strsplit(x[1], "/")[[1]]
  return(as.numeric(parts[1]))  
})

## Warning in FUN(X[[i]], ...): NAs introduced by coercion

Questionnaire$Arousal_4 <- sapply(strsplit(as.character(Questionnaire$C_4), ","), function(x) {
  parts <- strsplit(x[length(x)], "/")[[1]]
  
  if (length(parts) == 2) {
    return(as.numeric(parts[2]))  
  } else {
    return(NA)  
  }
})

Adjusting the Experiment Dataset

# Only include 'responsprocedure' trials & rename dataset
Experiment <- Experiment %>% filter(procedure == "responsprocedure")

# Renaming Colunmns
Experiment <- Experiment %>%
  rename(
    ID = subject,
    Block = session,
    Step = sub.trial.number,
    `Correct?` = feedback.ACC,
    Correct_Answer = feedback.CRESP,
    Given_Answer = feedback.RESP,
    Response_Time = feedback.RT
  )

# Adding a 'Trial' column (1-48)
Experiment <- Experiment %>%
  group_by(Block, ID) %>%
  mutate(Trial = rep(1:48, each = case_when(
    Block == 1  ~ 6,
    Block == 2  ~ 12,
    Block == 18 ~ 18,
    TRUE ~ n() / 48  # Default: dynamically calculate N per trial
  ))[1:n()]) %>%  # Ensures correct length
  ungroup()

## Warning: There were 90 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `Trial = ...[]`.
## ℹ In group 1: `Block = 1` `ID = 2`.
## Caused by warning in `rep()`:
## ! first element used of 'each' argument
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 89 remaining warnings.

# Adding a 'Trial_Accuracy' Column
Experiment <- Experiment %>% mutate(`Correct?` = as.numeric(`Correct?`), Response_Time = as.numeric(Response_Time))
Experiment <- Experiment %>%
  group_by(ID, Block, Trial) %>%
  mutate(Trial_Accuracy = sum(`Correct?`, na.rm = TRUE) / n()) %>%
  ungroup()

# Adding a 'Avg_Trial_RT' Column
Experiment <- Experiment %>%
  group_by(ID, Block, Trial) %>%
  mutate(Avg_Trial_RT = sum(Response_Time, na.rm = TRUE) / n()) %>%
  ungroup()

Investigating the Distribution of trial RT (before applying the filters)

# Check Individual Range of Trial RTs
range_per_participant_prior <- Experiment %>%
  group_by(ID) %>%
  summarise(
    Min_RT = min(Avg_Trial_RT, na.rm = TRUE),
    Max_RT = max(Avg_Trial_RT, na.rm = TRUE)
  )

# Check Distribution (Mean, Median, Global min max)
summary(Experiment$Avg_Trial_RT)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   68.92  367.67  502.46  575.65  716.77 5071.50

Filter Trials by RT

# Exclude trials where RT > 2 SDs above mean (within ID & Block))
Experiment_Clean <- Experiment %>%
  group_by(ID, Block) %>%
  mutate(RT_Mean = mean(Response_Time, na.rm = TRUE),
         RT_SD = sd(Response_Time, na.rm = TRUE)) %>%
  filter(Response_Time < RT_Mean + 2 * RT_SD) %>%  
  ungroup()


# Check Distribution (Mean, Median, Global min max)
summary(Experiment_Clean$Avg_Trial_RT)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   68.92  364.83  497.58  567.54  702.83 5071.50

Filter Trials by Accuracy

# Filter trials based on Accuracy (Remove if Accuracy < 0.8)
Experiment_Clean <- Experiment_Clean %>%filter(Trial_Accuracy >= 0.8) %>% mutate (ID=as.factor(ID))

# Check Distribution 
summary(Experiment_Clean$Trial_Accuracy)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.8333  1.0000  1.0000  0.9823  1.0000  1.0000

Comparing every IDs total trials before and after filtering

trial_count_per_ID_prior <- Experiment %>%
  group_by(ID) %>%
  summarise(Trial_Count = n())  

trial_count_per_ID_after <- Experiment_Clean %>%
  group_by(ID) %>%
  summarise(Trial_Count = n())  

print(trial_count_per_ID_prior)

## # A tibble: 19 × 2
##       ID Trial_Count
##    <int>       <int>
##  1     2        2880
##  2     3        2880
##  3     4        2880
##  4     5        2880
##  5     7        2880
##  6     8        2880
##  7    10        2880
##  8    11        2880
##  9    12         864
## 10    13        2880
## 11    14        2880
## 12    15        2880
## 13    16        2880
## 14    17        2880
## 15    18        2880
## 16    19        2880
## 17    20        2880
## 18    22        2880
## 19    23        2880

print(trial_count_per_ID_after)

## # A tibble: 19 × 2
##    ID    Trial_Count
##    <fct>       <int>
##  1 2            2437
##  2 3            2030
##  3 4            1964
##  4 5            1927
##  5 7            2202
##  6 8            2061
##  7 10           2045
##  8 11           1721
##  9 12            641
## 10 13           2266
## 11 14           2153
## 12 15           1907
## 13 16           2349
## 14 17           2559
## 15 18           1863
## 16 19           1311
## 17 20           1600
## 18 22           1522
## 19 23           2224

Merging the 2 Datasets

# Preparation
Questionnaire <- Questionnaire %>% mutate(ID = as.factor(ID),Sequence= as.factor(Sequence))

# Merging
Merged <- Experiment_Clean %>% mutate(ID = as.factor(ID)) 
Merged <- Merged %>% left_join(Questionnaire, by = "ID")
Merged <- Merged %>% mutate(Block= as.factor(Block), Step = as.factor(Step)) 

# Create new Variable BMI
Merged <- Merged %>%
  mutate(BMI = Weight / ((Height / 100) * (Height / 100)))

# Removing ID 12 (technical issues) and ID 19 (too inaccurate trials (<50% kept after filters))
Merged <- Merged %>% 
  filter(!ID %in% c(12, 19))

Descriptive Statistics Age

summary(Merged$Age)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   18.00   19.00   21.00   21.11   23.00   29.00

sd(Merged$Age)

## [1] 2.680631

Descriptive Statistics Footedness

absolute_frequencies_foot <- Merged %>%
  group_by(Footedness) %>%
  summarise(Count = n_distinct(ID))  
relative_frequencies_foot <- prop.table(absolute_frequencies_foot$Count)
absolute_frequencies_foot

## # A tibble: 2 × 2
##   Footedness Count
##        <int> <int>
## 1          1    16
## 2          2     1

relative_frequencies_foot

## [1] 0.94117647 0.05882353

Descriptive Statistics Gender

absolute_frequencies_gender <- Merged %>%
  group_by(Gender) %>%
  summarise(Count = n_distinct(ID))  
relative_frequencies_gender <- prop.table(absolute_frequencies_gender$Count)
absolute_frequencies_gender

## # A tibble: 2 × 2
##   Gender Count
##    <int> <int>
## 1      1     7
## 2      2    10

relative_frequencies_gender

## [1] 0.4117647 0.5882353

Descriptive Statistics Smoking

absolute_frequencies_smoking <- Merged %>%
  group_by(Smoking) %>%
  summarise(Count = n_distinct(ID))  
relative_frequencies_smoking <- prop.table(absolute_frequencies_smoking$Count)
absolute_frequencies_smoking

## # A tibble: 2 × 2
##   Smoking Count
##     <int> <int>
## 1       1     3
## 2       2    14

relative_frequencies_smoking

## [1] 0.1764706 0.8235294

Descriptive Statistics Alcohol (last 24 hours)

absolute_frequencies_alcohol <- Merged %>%
  group_by(Alcohol) %>%
  summarise(Count = n_distinct(ID))  
relative_frequencies_alcohol <- prop.table(absolute_frequencies_alcohol$Count)
absolute_frequencies_alcohol

## # A tibble: 2 × 2
##   Alcohol Count
##     <int> <int>
## 1       1     5
## 2       2    12

relative_frequencies_alcohol

## [1] 0.2941176 0.7058824

Reshaping the Dataset into Long Form

## Adding variable Avg_Block_RT to Merged
Merged <- Merged %>%
  mutate(Block = factor(Block, levels = c("B", "1", "2", "3", "4"))) %>%
  
group_by(ID, Block) %>%
mutate(Avg_Block_RT = mean(Avg_Trial_RT, na.rm = TRUE)) %>%
ungroup()


## Convert 'Block' to an ordered factor 
Merged <- Merged %>%
  mutate(Block = factor(Block, levels = c("B", "1", "2", "3", "4"), ordered = TRUE))

### Turning it long form
Merged_Long <- Merged %>%
  select(ID, 
         Q1_B, Q1_1, Q1_2, Q1_3, Q1_4,  # Mental Demand variables
         Q2_B, Q2_1, Q2_2, Q2_3, Q2_4,  # Physical Demand variables
         Q3_B, Q3_1, Q3_2, Q3_3, Q3_4,  # Temporal Demand variables
         Q4_B, Q4_1, Q4_2, Q4_3, Q4_4,  # Performance variables
         Q5_B, Q5_1, Q5_2, Q5_3, Q5_4,  # Effort variables
         Q6_B, Q6_1, Q6_2, Q6_3, Q6_4,  # Frustration variables
         Arousal_B, Arousal_1, Arousal_2, Arousal_3, Arousal_4,  # Arousal variables
         Pleasure_B, Pleasure_1, Pleasure_2, Pleasure_3, Pleasure_4) %>%
  
  # Convert the wide format to long format
  pivot_longer(cols = -c(ID), names_to = c("Measure", "Block_Num"), names_sep = "_", values_to = "Value") %>%
  
  # Set block as a factor, keeping the correct block order
  mutate(Block = factor(Block_Num, levels = c("B", "1", "2", "3", "4"), ordered = TRUE)) %>%
  
  # Remove the Block_Num column (because Block is now the correct one)
  select(-Block_Num) %>%
  
  # Group by ID and Block, and calculate the mean Value for each
  group_by(ID, Block, Measure) %>%
  summarise(Value = mean(Value, na.rm = TRUE), .groups = "drop") %>%
  
  # Pivot back to wide format with one row per ID and Block
  pivot_wider(names_from = Measure, values_from = Value) %>%
  
  # Join Avg_Block_RT from the Merged dataset
  left_join(Merged %>%
              select(ID, Block, Avg_Block_RT) %>%
              distinct(), 
            by = c("ID", "Block"))

## Rename Columns
Merged_Long <- Merged_Long %>%
  rename(
    Mental_Demand = Q1,
    Physical_Demand = Q2,
    Temporal_Demand = Q3,
    Performance = Q4,
    Effort = Q5,
    Frustration = Q6
  )

Plotting & Analysis

Mental Demand by Block

LMER_MD <- lmer(Mental_Demand ~ Block  + (1 | ID) , data = Merged_Long)

Anova(LMER_MD)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Mental_Demand
##        Chisq Df Pr(>Chisq)    
## Block 39.617  4  5.193e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Posthoc: emmeans

EMM_MD <- emmeans(LMER_MD, pairwise ~ Block)
print(EMM_MD)

## $emmeans
##  Block emmean   SE   df lower.CL upper.CL
##  B       3.65 1.02 61.4     1.61     5.68
##  1       6.71 1.02 61.4     4.67     8.74
##  2       9.24 1.02 61.4     7.20    11.27
##  3      10.53 1.02 61.4     8.50    12.56
##  4       9.27 1.04 63.6     7.18    11.35
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast estimate   SE   df t.ratio p.value
##  B - 1     -3.0588 1.23 63.0  -2.488  0.1064
##  B - 2     -5.5882 1.23 63.0  -4.545  0.0002
##  B - 3     -6.8824 1.23 63.0  -5.597  <.0001
##  B - 4     -5.6205 1.25 63.4  -4.489  0.0003
##  1 - 2     -2.5294 1.23 63.0  -2.057  0.2516
##  1 - 3     -3.8235 1.23 63.0  -3.110  0.0227
##  1 - 4     -2.5616 1.25 63.4  -2.046  0.2565
##  2 - 3     -1.2941 1.23 63.0  -1.053  0.8297
##  2 - 4     -0.0322 1.25 63.4  -0.026  1.0000
##  3 - 4      1.2619 1.25 63.4   1.008  0.8508
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 5 estimates

Plot

#Effects
ae.m.MD <- allEffects(LMER_MD)
ae.m.MD.df <- as.data.frame(ae.m.MD[[1]])

#Confidence Intervals
#change conf interval to 83%, match p = .05
ae.m.MD.df$l83 <- ae.m.MD.df$fit - 1.3722 * ae.m.MD.df$se
ae.m.MD.df$u83 <- ae.m.MD.df$fit + 1.3722 * ae.m.MD.df$se

## Effects Plot
plot.MD <- ggplot(ae.m.MD.df, aes(x=Block, y=fit, ymin=l83, ymax=u83))+
  geom_point(aes(color = Block), size=3) +
  geom_errorbar(aes(ymin=l83, ymax=u83, color= Block), width=.1, size=1) +
  geom_line() +
  ylab("Mental Demand")+
  xlab("Block")+
  ggtitle("Mental Demand across Blocks")+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

print(plot.MD)

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

# Plotting Residuals and Assumption of Normality

residuals <- resid(LMER_MD)
fitted <- fitted(LMER_MD)

par(mfrow = c(1, 2))
hist(residuals, main = "Histogram of Residuals", xlab = "Residuals", breaks = 20)
qqnorm(residuals)
qqline(residuals, col = "red")

summary(LMER_MD, ddf="Satterthwaite")

## Warning in summary.merMod(LMER_MD, ddf = "Satterthwaite"): additional arguments
## ignored

## Linear mixed model fit by REML ['lmerMod']
## Formula: Mental_Demand ~ Block + (1 | ID)
##    Data: Merged_Long
## 
## REML criterion at convergence: 458.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8081 -0.6341 -0.1137  0.6339  2.5240 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept)  4.732   2.175   
##  Residual             12.850   3.585   
## Number of obs: 84, groups:  ID, 17
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  7.87703    0.65703  11.989
## Block.L      4.76380    0.88182   5.402
## Block.Q     -2.63967    0.87829  -3.005
## Block.C     -0.64087    0.87253  -0.734
## Block^4     -0.07349    0.86986  -0.084
## 
## Correlation of Fixed Effects:
##         (Intr) Blck.L Blck.Q Blck.C
## Block.L 0.012                      
## Block.Q 0.010  0.024               
## Block.C 0.006  0.014  0.012        
## Block^4 0.002  0.005  0.005  0.003

Physical Demand by Block

LMER_PD <- lmer(Physical_Demand ~ Block  + (1 | ID) , data = Merged_Long)

Anova(LMER_PD)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Physical_Demand
##        Chisq Df Pr(>Chisq)    
## Block 28.628  4  9.303e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#Posthoc: emmeans

EMM_PD <- emmeans(LMER_PD, pairwise ~ Block)
print(EMM_PD)

## $emmeans
##  Block emmean    SE   df lower.CL upper.CL
##  B       3.29 0.822 44.3     1.64     4.95
##  1       4.41 0.822 44.3     2.75     6.07
##  2       6.47 0.822 44.3     4.81     8.13
##  3       7.35 0.822 44.3     5.70     9.01
##  4       6.00 0.822 44.3     4.34     7.66
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast estimate    SE df t.ratio p.value
##  B - 1      -1.118 0.864 64  -1.294  0.6957
##  B - 2      -3.176 0.864 64  -3.678  0.0043
##  B - 3      -4.059 0.864 64  -4.700  0.0001
##  B - 4      -2.706 0.864 64  -3.133  0.0212
##  1 - 2      -2.059 0.864 64  -2.384  0.1330
##  1 - 3      -2.941 0.864 64  -3.406  0.0097
##  1 - 4      -1.588 0.864 64  -1.839  0.3607
##  2 - 3      -0.882 0.864 64  -1.022  0.8444
##  2 - 4       0.471 0.864 64   0.545  0.9822
##  3 - 4       1.353 0.864 64   1.567  0.5239
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 5 estimates

Plot

#Effects
ae.m.PD <- allEffects(LMER_PD)
ae.m.PD.df <- as.data.frame(ae.m.PD[[1]])

#Confidence Intervals
#change conf interval to 83%, match p = .05
ae.m.PD.df$l83 <- ae.m.PD.df$fit - 1.3722 * ae.m.PD.df$se
ae.m.PD.df$u83 <- ae.m.PD.df$fit + 1.3722 * ae.m.PD.df$se

## Effects Plot
plot.PD <- ggplot(ae.m.PD.df, aes(x=Block, y=fit, ymin=l83, ymax=u83))+
  geom_point(aes(color = Block), size=3) +
  geom_errorbar(aes(ymin=l83, ymax=u83, color= Block), width=.1, size=1) +
  geom_line() +
  ylab("Physical Demand")+
  xlab("Block")+
  ggtitle("Physical Demand across Blocks")+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))

print(plot.PD)

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

# Plotting Residuals and Assumption of Normality

residuals <- resid(LMER_PD)
fitted <- fitted(LMER_PD)

par(mfrow = c(1, 2))
hist(residuals, main = "Histogram of Residuals", xlab = "Residuals", breaks = 20)
qqnorm(residuals)
qqline(residuals, col = "red")

summary(LMER_PD, ddf="Satterthwaite")

## Warning in summary.merMod(LMER_PD, ddf = "Satterthwaite"): additional arguments
## ignored

## Linear mixed model fit by REML ['lmerMod']
## Formula: Physical_Demand ~ Block + (1 | ID)
##    Data: Merged_Long
## 
## REML criterion at convergence: 416.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6425 -0.6924 -0.0573  0.4966  3.3002 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 5.158    2.271   
##  Residual             6.339    2.518   
## Number of obs: 85, groups:  ID, 17
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   5.5059     0.6148   8.956
## Block.L       2.6414     0.6107   4.326
## Block.Q      -1.6350     0.6107  -2.677
## Block.C      -1.0045     0.6107  -1.645
## Block^4       0.1266     0.6107   0.207
## 
## Correlation of Fixed Effects:
##         (Intr) Blck.L Blck.Q Blck.C
## Block.L 0.000                      
## Block.Q 0.000  0.000               
## Block.C 0.000  0.000  0.000        
## Block^4 0.000  0.000  0.000  0.000

Temporal Demand by Block

LMER_TD <- lmer(Temporal_Demand ~ Block  + (1 | ID) , data = Merged_Long)

Anova(LMER_TD)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Temporal_Demand
##        Chisq Df Pr(>Chisq)   
## Block 15.992  4    0.00303 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#Posthoc: emmeans

EMM_TD <- emmeans(LMER_TD, pairwise ~ Block)
print(EMM_TD)

## $emmeans
##  Block emmean    SE   df lower.CL upper.CL
##  B       3.41 0.781 49.4     1.84     4.98
##  1       5.53 0.781 49.4     3.96     7.10
##  2       6.41 0.781 49.4     4.84     7.98
##  3       6.35 0.781 49.4     4.78     7.92
##  4       5.29 0.781 49.4     3.72     6.86
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast estimate   SE df t.ratio p.value
##  B - 1     -2.1176 0.86 64  -2.463  0.1123
##  B - 2     -3.0000 0.86 64  -3.489  0.0076
##  B - 3     -2.9412 0.86 64  -3.420  0.0093
##  B - 4     -1.8824 0.86 64  -2.189  0.1971
##  1 - 2     -0.8824 0.86 64  -1.026  0.8424
##  1 - 3     -0.8235 0.86 64  -0.958  0.8729
##  1 - 4      0.2353 0.86 64   0.274  0.9987
##  2 - 3      0.0588 0.86 64   0.068  1.0000
##  2 - 4      1.1176 0.86 64   1.300  0.6923
##  3 - 4      1.0588 0.86 64   1.231  0.7333
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 5 estimates

Plot

#Effects
ae.m.TD <- allEffects(LMER_TD)
ae.m.TD.df <- as.data.frame(ae.m.TD[[1]])

#Confidence Intervals
#change conf interval to 83%, match p = .05
ae.m.TD.df$l83 <- ae.m.TD.df$fit - 1.3722 * ae.m.TD.df$se
ae.m.TD.df$u83 <- ae.m.TD.df$fit + 1.3722 * ae.m.TD.df$se

## Effects Plot
plot.TD <- ggplot(ae.m.TD.df, aes(x=Block, y=fit, ymin=l83, ymax=u83))+
  geom_point(aes(color = Block), size=3) +
  geom_errorbar(aes(ymin=l83, ymax=u83, color= Block), width=.1, size=1) +
  geom_line() +
  ylab("Temporal Demand")+
  xlab("Block")+
  ggtitle("Temporal Demand across Blocks")+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))

print(plot.TD)

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

# Plotting Residuals and Assumption of Normality

residuals <- resid(LMER_TD)
fitted <- fitted(LMER_TD)

par(mfrow = c(1, 2))
hist(residuals, main = "Histogram of Residuals", xlab = "Residuals", breaks = 20)
qqnorm(residuals)
qqline(residuals, col = "red")

summary(LMER_TD, ddf="Satterthwaite")

## Warning in summary.merMod(LMER_TD, ddf = "Satterthwaite"): additional arguments
## ignored

## Linear mixed model fit by REML ['lmerMod']
## Formula: Temporal_Demand ~ Block + (1 | ID)
##    Data: Merged_Long
## 
## REML criterion at convergence: 413
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.7316 -0.5818 -0.1126  0.4646  3.0171 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 4.088    2.022   
##  Residual             6.286    2.507   
## Number of obs: 85, groups:  ID, 17
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  5.40000    0.56072   9.630
## Block.L      1.45093    0.60807   2.386
## Block.Q     -1.94943    0.60807  -3.206
## Block.C      0.07441    0.60807   0.122
## Block^4     -0.04218    0.60807  -0.069
## 
## Correlation of Fixed Effects:
##         (Intr) Blck.L Blck.Q Blck.C
## Block.L 0.000                      
## Block.Q 0.000  0.000               
## Block.C 0.000  0.000  0.000        
## Block^4 0.000  0.000  0.000  0.000

Performance by Block

LMER_P <- lmer(Performance ~ Block  + (1 | ID) , data = Merged_Long)

Anova(LMER_P)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Performance
##        Chisq Df Pr(>Chisq)    
## Block 24.965  4  5.114e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Posthoc: emmeans

EMM_P <- emmeans(LMER_P, pairwise ~ Block)
print(EMM_P)

## $emmeans
##  Block emmean  SE   df lower.CL upper.CL
##  B       6.65 1.2 67.2     4.24     9.05
##  1      10.88 1.2 67.2     8.48    13.29
##  2      12.82 1.2 67.2    10.42    15.23
##  3      11.88 1.2 67.2     9.48    14.29
##  4      13.29 1.2 67.2    10.89    15.70
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast estimate  SE df t.ratio p.value
##  B - 1      -4.235 1.5 64  -2.814  0.0491
##  B - 2      -6.176 1.5 64  -4.104  0.0011
##  B - 3      -5.235 1.5 64  -3.479  0.0078
##  B - 4      -6.647 1.5 64  -4.417  0.0004
##  1 - 2      -1.941 1.5 64  -1.290  0.6983
##  1 - 3      -1.000 1.5 64  -0.664  0.9632
##  1 - 4      -2.412 1.5 64  -1.603  0.5013
##  2 - 3       0.941 1.5 64   0.625  0.9705
##  2 - 4      -0.471 1.5 64  -0.313  0.9979
##  3 - 4      -1.412 1.5 64  -0.938  0.8810
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 5 estimates

Plot

#Effects
ae.m.P <- allEffects(LMER_P)
ae.m.P.df <- as.data.frame(ae.m.P[[1]])

#Confidence Intervals
#change conf interval to 83%, match p = .05
ae.m.P.df$l83 <- ae.m.P.df$fit - 1.3722 * ae.m.P.df$se
ae.m.P.df$u83 <- ae.m.P.df$fit + 1.3722 * ae.m.P.df$se

## Effects Plot
plot.P <- ggplot(ae.m.P.df, aes(x=Block, y=fit, ymin=l83, ymax=u83))+
  geom_point(aes(color = Block), size=3) +
  geom_errorbar(aes(ymin=l83, ymax=u83, color= Block), width=.1, size=1) +
  geom_line() +
  ylab("Performance")+
  xlab("Block")+
  ggtitle("Performance across Blocks")+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))

print(plot.P)

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

# Plotting Residuals and Assumption of Normality

residuals <- resid(LMER_P)
fitted <- fitted(LMER_P)

par(mfrow = c(1, 2))
hist(residuals, main = "Histogram of Residuals", xlab = "Residuals", breaks = 20)
qqnorm(residuals)
qqline(residuals, col = "red")

summary(LMER_P, ddf="Satterthwaite")

## Warning in summary.merMod(LMER_P, ddf = "Satterthwaite"): additional arguments
## ignored

## Linear mixed model fit by REML ['lmerMod']
## Formula: Performance ~ Block + (1 | ID)
##    Data: Merged_Long
## 
## REML criterion at convergence: 493.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.9376 -0.6808 -0.1318  0.6269  3.0193 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept)  5.39    2.322   
##  Residual             19.25    4.388   
## Number of obs: 85, groups:  ID, 17
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  11.1059     0.7373  15.064
## Block.L       4.5202     1.0642   4.248
## Block.Q      -2.2796     1.0642  -2.142
## Block.C       1.4695     1.0642   1.381
## Block^4       0.6960     1.0642   0.654
## 
## Correlation of Fixed Effects:
##         (Intr) Blck.L Blck.Q Blck.C
## Block.L 0.000                      
## Block.Q 0.000  0.000               
## Block.C 0.000  0.000  0.000        
## Block^4 0.000  0.000  0.000  0.000

Effort across Blocks

LMER_E <- lmer(Effort ~ Block  + (1 | ID) , data = Merged_Long)

Anova(LMER_E)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Effort
##        Chisq Df Pr(>Chisq)    
## Block 36.114  4  2.741e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Posthoc: emmeans

EMM_E <- emmeans(LMER_E, pairwise ~ Block)
print(EMM_E)

## $emmeans
##  Block emmean SE df lower.CL upper.CL
##  B       3.82  1 69     1.83     5.82
##  1       7.76  1 69     5.77     9.76
##  2      10.12  1 69     8.12    12.11
##  3      10.53  1 69     8.54    12.52
##  4       8.88  1 69     6.89    10.88
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast estimate   SE df t.ratio p.value
##  B - 1      -3.941 1.27 64  -3.115  0.0223
##  B - 2      -6.294 1.27 64  -4.974  <.0001
##  B - 3      -6.706 1.27 64  -5.300  <.0001
##  B - 4      -5.059 1.27 64  -3.998  0.0015
##  1 - 2      -2.353 1.27 64  -1.860  0.3494
##  1 - 3      -2.765 1.27 64  -2.185  0.1986
##  1 - 4      -1.118 1.27 64  -0.883  0.9020
##  2 - 3      -0.412 1.27 64  -0.325  0.9975
##  2 - 4       1.235 1.27 64   0.976  0.8649
##  3 - 4       1.647 1.27 64   1.302  0.6911
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 5 estimates

Plot

#Effects
ae.m.E <- allEffects(LMER_E)
ae.m.E.df <- as.data.frame(ae.m.E[[1]])

#Confidence Intervals
#change conf interval to 83%, match p = .05
ae.m.E.df$l83 <- ae.m.E.df$fit - 1.3722 * ae.m.E.df$se
ae.m.E.df$u83 <- ae.m.E.df$fit + 1.3722 * ae.m.E.df$se

## Effects Plot
plot.E <- ggplot(ae.m.E.df, aes(x=Block, y=fit, ymin=l83, ymax=u83))+
  geom_point(aes(color = Block), size=3) +
  geom_errorbar(aes(ymin=l83, ymax=u83, color= Block), width=.1, size=1) +
  geom_line() +
  ylab("Effort")+
  xlab("Block")+
  ggtitle("Effort across Blocks")+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))

print(plot.E)

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

# Plotting Residuals and Assumption of Normality

residuals <- resid(LMER_E)
fitted <- fitted(LMER_E)

par(mfrow = c(1, 2))
hist(residuals, main = "Histogram of Residuals", xlab = "Residuals", breaks = 20)
qqnorm(residuals)
qqline(residuals, col = "red")

summary(LMER_E, ddf="Satterthwaite")

## Warning in summary.merMod(LMER_E, ddf = "Satterthwaite"): additional arguments
## ignored

## Linear mixed model fit by REML ['lmerMod']
## Formula: Effort ~ Block + (1 | ID)
##    Data: Merged_Long
## 
## REML criterion at convergence: 464.6
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.83578 -0.68975  0.00585  0.42436  2.16060 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept)  3.383   1.839   
##  Residual             13.608   3.689   
## Number of obs: 85, groups:  ID, 17
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  8.22353    0.59924  13.723
## Block.L      4.07376    0.89471   4.553
## Block.Q     -3.50584    0.89471  -3.918
## Block.C     -0.14881    0.89471  -0.166
## Block^4      0.02812    0.89471   0.031
## 
## Correlation of Fixed Effects:
##         (Intr) Blck.L Blck.Q Blck.C
## Block.L 0.000                      
## Block.Q 0.000  0.000               
## Block.C 0.000  0.000  0.000        
## Block^4 0.000  0.000  0.000  0.000

Frustration by Block

LMER_F <- lmer(Frustration ~ Block  + (1 | ID) , data = Merged_Long)

Anova(LMER_F)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Frustration
##        Chisq Df Pr(>Chisq)   
## Block 16.352  4   0.002581 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Posthoc: emmeans

EMM_F <- emmeans(LMER_F, pairwise ~ Block)
print(EMM_F)

## $emmeans
##  Block emmean    SE   df lower.CL upper.CL
##  B       4.24 0.993 77.2     2.26     6.21
##  1       6.76 0.993 77.2     4.79     8.74
##  2       7.47 0.993 77.2     5.49     9.45
##  3       9.53 0.993 77.2     7.55    11.51
##  4       6.41 0.993 77.2     4.43     8.39
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast estimate   SE df t.ratio p.value
##  B - 1      -2.529 1.34 64  -1.893  0.3315
##  B - 2      -3.235 1.34 64  -2.422  0.1228
##  B - 3      -5.294 1.34 64  -3.962  0.0017
##  B - 4      -2.176 1.34 64  -1.629  0.4848
##  1 - 2      -0.706 1.34 64  -0.528  0.9841
##  1 - 3      -2.765 1.34 64  -2.069  0.2460
##  1 - 4       0.353 1.34 64   0.264  0.9989
##  2 - 3      -2.059 1.34 64  -1.541  0.5402
##  2 - 4       1.059 1.34 64   0.792  0.9319
##  3 - 4       3.118 1.34 64   2.333  0.1478
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 5 estimates

Plot

#Effects
ae.m.F <- allEffects(LMER_F)
ae.m.F.df <- as.data.frame(ae.m.F[[1]])

#Confidence Intervals
#change conf interval to 83%, match p = .05
ae.m.F.df$l83 <- ae.m.F.df$fit - 1.3722 * ae.m.F.df$se
ae.m.F.df$u83 <- ae.m.F.df$fit + 1.3722 * ae.m.F.df$se

## Effects Plot
plot.F <- ggplot(ae.m.F.df, aes(x=Block, y=fit, ymin=l83, ymax=u83))+
  geom_point(aes(color = Block), size=3) +
  geom_errorbar(aes(ymin=l83, ymax=u83, color= Block), width=.1, size=1) +
  geom_line() +
  ylab("Frustration")+
  xlab("Block")+
  ggtitle("Frustration across Blocks")+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))

print(plot.F)

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

# Plotting Residuals and Assumption of Normality

residuals <- resid(LMER_F)
fitted <- fitted(LMER_F)

par(mfrow = c(1, 2))
hist(residuals, main = "Histogram of Residuals", xlab = "Residuals", breaks = 20)
qqnorm(residuals)
qqline(residuals, col = "red")

summary(LMER_F, ddf="Satterthwaite")

## Warning in summary.merMod(LMER_F, ddf = "Satterthwaite"): additional arguments
## ignored

## Linear mixed model fit by REML ['lmerMod']
## Formula: Frustration ~ Block + (1 | ID)
##    Data: Merged_Long
## 
## REML criterion at convergence: 467.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.5633 -0.7176 -0.1689  0.5560  2.5543 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept)  1.586   1.259   
##  Residual             15.173   3.895   
## Number of obs: 85, groups:  ID, 17
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   6.8824     0.5213  13.202
## Block.L       2.2508     0.9447   2.382
## Block.Q      -2.6569     0.9447  -2.812
## Block.C      -1.0603     0.9447  -1.122
## Block^4      -1.1601     0.9447  -1.228
## 
## Correlation of Fixed Effects:
##         (Intr) Blck.L Blck.Q Blck.C
## Block.L 0.000                      
## Block.Q 0.000  0.000               
## Block.C 0.000  0.000  0.000        
## Block^4 0.000  0.000  0.000  0.000

Pleasure by Block

LMER_p <- lmer(Pleasure ~ Block  + (1 | ID) , data = Merged_Long)

Anova(LMER_p)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Pleasure
##        Chisq Df Pr(>Chisq)
## Block 1.8207  4     0.7687

Plot

#Effects
ae.m.p <- allEffects(LMER_p)
ae.m.p.df <- as.data.frame(ae.m.p[[1]])

#Confidence Intervals
#change conf interval to 83%, match p = .05
ae.m.p.df$l83 <- ae.m.p.df$fit - 1.3722 * ae.m.p.df$se
ae.m.p.df$u83 <- ae.m.p.df$fit + 1.3722 * ae.m.p.df$se

## Effects Plot
plot.p <- ggplot(ae.m.p.df, aes(x=Block, y=fit, ymin=l83, ymax=u83))+
  geom_point(aes(color = Block), size=3) +
  geom_errorbar(aes(ymin=l83, ymax=u83, color= Block), width=.1, size=1) +
  geom_line() +
  ylab("Pleasure")+
  xlab("Block")+
  ggtitle("Pleasure across Blocks")+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))

print(plot.p)

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

# Plotting Residuals and Assumption of Normality

residuals <- resid(LMER_p)
fitted <- fitted(LMER_p)

par(mfrow = c(1, 2))
hist(residuals, main = "Histogram of Residuals", xlab = "Residuals", breaks = 20)
qqnorm(residuals)
qqline(residuals, col = "red")

summary(LMER_p, ddf="Satterthwaite")

## Warning in summary.merMod(LMER_p, ddf = "Satterthwaite"): additional arguments
## ignored

## Linear mixed model fit by REML ['lmerMod']
## Formula: Pleasure ~ Block + (1 | ID)
##    Data: Merged_Long
## 
## REML criterion at convergence: 317.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.5394 -0.5048  0.1896  0.6025  1.6368 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 0.7916   0.8897  
##  Residual             2.1542   1.4677  
## Number of obs: 84, groups:  ID, 17
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   6.1122     0.2688  22.736
## Block.L       0.2989     0.3611   0.828
## Block.Q       0.3470     0.3596   0.965
## Block.C       0.1030     0.3573   0.288
## Block^4      -0.1439     0.3562  -0.404
## 
## Correlation of Fixed Effects:
##         (Intr) Blck.L Blck.Q Blck.C
## Block.L 0.012                      
## Block.Q 0.010  0.024               
## Block.C 0.006  0.014  0.012        
## Block^4 0.002  0.005  0.005  0.003

Arousal by Block

LMER_a <- lmer(Arousal ~ Block  + (1 | ID) , data = Merged_Long)

Anova(LMER_a)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Arousal
##        Chisq Df Pr(>Chisq)  
## Block 13.087  4    0.01086 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Posthoc: emmeans

EMM_a <- emmeans(LMER_a, pairwise ~ Block)
print(EMM_a)

## $emmeans
##  Block emmean    SE   df lower.CL upper.CL
##  B       5.41 0.432 46.2     4.54     6.28
##  1       6.06 0.432 46.2     5.19     6.93
##  2       6.35 0.432 46.2     5.48     7.22
##  3       5.35 0.432 46.2     4.48     6.22
##  4       4.84 0.442 48.5     3.95     5.73
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast estimate    SE   df t.ratio p.value
##  B - 1     -0.6471 0.464 63.0  -1.395  0.6334
##  B - 2     -0.9412 0.464 63.0  -2.028  0.2647
##  B - 3      0.0588 0.464 63.0   0.127  0.9999
##  B - 4      0.5705 0.473 63.2   1.207  0.7473
##  1 - 2     -0.2941 0.464 63.0  -0.634  0.9690
##  1 - 3      0.7059 0.464 63.0   1.521  0.5528
##  1 - 4      1.2176 0.473 63.2   2.576  0.0872
##  2 - 3      1.0000 0.464 63.0   2.155  0.2104
##  2 - 4      1.5117 0.473 63.2   3.198  0.0178
##  3 - 4      0.5117 0.473 63.2   1.083  0.8148
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 5 estimates

Plot

#Effects
ae.m.a <- allEffects(LMER_a)
ae.m.a.df <- as.data.frame(ae.m.a[[1]])

#Confidence Intervals
#change conf interval to 83%, match p = .05
ae.m.a.df$l83 <- ae.m.a.df$fit - 1.3722 * ae.m.a.df$se
ae.m.a.df$u83 <- ae.m.a.df$fit + 1.3722 * ae.m.a.df$se

## Effects Plot
plot.a <- ggplot(ae.m.a.df, aes(x=Block, y=fit, ymin=l83, ymax=u83))+
  geom_point(aes(color = Block), size=3) +
  geom_errorbar(aes(ymin=l83, ymax=u83, color= Block), width=.1, size=1) +
  geom_line() +
  ylab("Arousal")+
  xlab("Block")+
  ggtitle("Arousal across Blocks")+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))

print(plot.a)

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

# Plotting Residuals and Assumption of Normality

residuals <- resid(LMER_a)
fitted <- fitted(LMER_a)

par(mfrow = c(1, 2))
hist(residuals, main = "Histogram of Residuals", xlab = "Residuals", breaks = 20)
qqnorm(residuals)
qqline(residuals, col = "red")

summary(LMER_a, ddf="Satterthwaite")

## Warning in summary.merMod(LMER_a, ddf = "Satterthwaite"): additional arguments
## ignored

## Linear mixed model fit by REML ['lmerMod']
## Formula: Arousal ~ Block + (1 | ID)
##    Data: Merged_Long
## 
## REML criterion at convergence: 312.2
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.12956 -0.53294  0.00254  0.79522  2.10447 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 1.348    1.161   
##  Residual             1.830    1.353   
## Number of obs: 84, groups:  ID, 17
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   5.6035     0.3180  17.618
## Block.L      -0.5840     0.3329  -1.754
## Block.Q      -0.9652     0.3316  -2.911
## Block.C       0.2660     0.3293   0.808
## Block^4       0.3255     0.3283   0.992
## 
## Correlation of Fixed Effects:
##         (Intr) Blck.L Blck.Q Blck.C
## Block.L 0.010                      
## Block.Q 0.008  0.024               
## Block.C 0.005  0.015  0.012        
## Block^4 0.002  0.006  0.005  0.003

RT by Block

LMER_RT <- lmer(Response_Time ~ Block  + (1 | ID) , data = Merged)

Anova(LMER_RT)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Response_Time
##        Chisq Df Pr(>Chisq)    
## Block 262.96  3  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Posthoc: emmeans

EMM_RT <- emmeans(LMER_RT, pairwise ~ Block)

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

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

print(EMM_RT)

## $emmeans
##  Block emmean   SE  df asymp.LCL asymp.UCL
##  1        480 50.6 Inf       381       579
##  2        466 50.5 Inf       367       565
##  3        487 50.5 Inf       388       586
##  4        436 50.5 Inf       337       535
## 
## Degrees-of-freedom method: asymptotic 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast        estimate   SE  df z.ratio p.value
##  Block1 - Block2     14.4 4.18 Inf   3.449  0.0032
##  Block1 - Block3     -7.1 4.00 Inf  -1.776  0.2851
##  Block1 - Block4     43.9 4.13 Inf  10.624  <.0001
##  Block2 - Block3    -21.5 3.31 Inf  -6.496  <.0001
##  Block2 - Block4     29.4 3.47 Inf   8.473  <.0001
##  Block3 - Block4     51.0 3.26 Inf  15.632  <.0001
## 
## Degrees-of-freedom method: asymptotic 
## P value adjustment: tukey method for comparing a family of 4 estimates

Plot

#Effects
ae.m.RT <- allEffects(LMER_RT)
ae.m.RT.df <- as.data.frame(ae.m.RT[[1]])

#Confidence Intervals
#change conf interval to 83%, match p = .05
ae.m.RT.df$l83 <- ae.m.RT.df$fit - 1.3722 * ae.m.RT.df$se
ae.m.RT.df$u83 <- ae.m.RT.df$fit + 1.3722 * ae.m.RT.df$se

## Effects Plot
plot.RT <- ggplot(ae.m.RT.df, aes(x=Block, y=fit, ymin=l83, ymax=u83))+
  geom_point(aes(color = Block), size=3) +
  geom_errorbar(aes(ymin=l83, ymax=u83, color= Block), width=.1, size=1) +
  geom_line() +
  ylab("Reaction Time (in ms)")+
  xlab("Block")+
  ggtitle("Reaction Time across Blocks")+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))

print(plot.RT)

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

# Plotting Residuals and Assumption of Normality

residuals <- resid(LMER_RT)
fitted <- fitted(LMER_RT)

par(mfrow = c(1, 2))
hist(residuals, main = "Histogram of Residuals", xlab = "Residuals", breaks = 20)
qqnorm(residuals)
qqline(residuals, col = "red")

summary(LMER_RT, ddf="Satterthwaite")

## Warning in summary.merMod(LMER_RT, ddf = "Satterthwaite"): additional arguments
## ignored

## Linear mixed model fit by REML ['lmerMod']
## Formula: Response_Time ~ Block + (1 | ID)
##    Data: Merged
## 
## REML criterion at convergence: 396519.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.4680 -0.6003 -0.2264  0.3613 13.0801 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 43257    208.0   
##  Residual             45866    214.2   
## Number of obs: 29210, groups:  ID, 17
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  467.173     50.461   9.258
## Block.L      -24.608      2.866  -8.586
## Block.Q      -18.273      2.654  -6.886
## Block.C      -24.244      2.407 -10.072
## 
## Correlation of Fixed Effects:
##         (Intr) Blck.L Blck.Q
## Block.L -0.007              
## Block.Q  0.005 -0.216       
## Block.C  0.000  0.132 -0.171

RT by Mental Demand

LMER_RT_MD <- lmer(Avg_Block_RT ~ Mental_Demand  + (1 | ID) , data = Merged_Long)

summary(LMER_RT_MD)

## Linear mixed model fit by REML ['lmerMod']
## Formula: Avg_Block_RT ~ Mental_Demand + (1 | ID)
##    Data: Merged_Long
## 
## REML criterion at convergence: 827.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.5890 -0.4846 -0.1370  0.5155  2.5391 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 46925    216.62  
##  Residual              7491     86.55  
## Number of obs: 67, groups:  ID, 17
## 
## Fixed effects:
##               Estimate Std. Error t value
## (Intercept)    493.082     60.199   8.191
## Mental_Demand    1.140      3.066   0.372
## 
## Correlation of Fixed Effects:
##             (Intr)
## Mental_Dmnd -0.455

Plot

ggplot(Merged_Long, aes(x = Mental_Demand, y = Avg_Block_RT)) +
  geom_point(color = "blue", size = 3, alpha = 0.7) +  
  geom_smooth(method = "lm", color = "red", size = 1.2, linetype = "solid") +  
  labs(title = "Reaction Time by Mental Demand",
       x = "Mental Demand",
       y = "Average Block Reaction Time (ms)") +
  theme_minimal() +
  theme(legend.position = "none") 

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 18 rows containing non-finite outside the scale range
## (`stat_smooth()`).

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

# Plotting Residuals and Assumption of Normality

residuals <- resid(LMER_RT_MD)
fitted <- fitted(LMER_RT_MD)

par(mfrow = c(1, 2))
hist(residuals, main = "Histogram of Residuals", xlab = "Residuals", breaks = 20)
qqnorm(residuals)
qqline(residuals, col = "red")

## RT by Physical Demand

LMER_RT_PD <- lmer(Avg_Block_RT ~ Physical_Demand  + (1 | ID) , data = Merged_Long)

summary(LMER_RT_PD)

## Linear mixed model fit by REML ['lmerMod']
## Formula: Avg_Block_RT ~ Physical_Demand + (1 | ID)
##    Data: Merged_Long
## 
## REML criterion at convergence: 838.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.3834 -0.5229 -0.1671  0.4941  2.2704 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 45221    212.65  
##  Residual              7651     87.47  
## Number of obs: 68, groups:  ID, 17
## 
## Fixed effects:
##                 Estimate Std. Error t value
## (Intercept)      530.572     61.061   8.689
## Physical_Demand   -4.869      5.103  -0.954
## 
## Correlation of Fixed Effects:
##             (Intr)
## Physcl_Dmnd -0.506

Plot

ggplot(Merged_Long, aes(x = Physical_Demand, y = Avg_Block_RT)) +
  geom_point(color = "blue", size = 3, alpha = 0.7) +  
  geom_smooth(method = "lm", color = "red", size = 1.2, linetype = "solid") +  
  labs(title = "Reaction Time by Physical Demand",
       x = "Physical Demand",
       y = "Average Block Reaction Time (ms)") +
  theme_minimal() +
  theme(legend.position = "none") 

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 17 rows containing non-finite outside the scale range
## (`stat_smooth()`).

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

# Plotting Residuals and Assumption of Normality

residuals <- resid(LMER_RT_PD)
fitted <- fitted(LMER_RT_PD)

par(mfrow = c(1, 2))
hist(residuals, main = "Histogram of Residuals", xlab = "Residuals", breaks = 20)
qqnorm(residuals)
qqline(residuals, col = "red")

## RT by Temporal Demand

LMER_RT_TD <- lmer(Avg_Block_RT ~ Temporal_Demand  + (1 | ID) , data = Merged_Long)

summary(LMER_RT_TD)

## Linear mixed model fit by REML ['lmerMod']
## Formula: Avg_Block_RT ~ Temporal_Demand + (1 | ID)
##    Data: Merged_Long
## 
## REML criterion at convergence: 836
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.23051 -0.56096 -0.07389  0.55737  2.29452 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 44480    210.90  
##  Residual              7318     85.54  
## Number of obs: 68, groups:  ID, 17
## 
## Fixed effects:
##                 Estimate Std. Error t value
## (Intercept)      559.954     60.945   9.188
## Temporal_Demand   -9.985      5.336  -1.871
## 
## Correlation of Fixed Effects:
##             (Intr)
## Temprl_Dmnd -0.516

Plot

ggplot(Merged_Long, aes(x = Temporal_Demand, y = Avg_Block_RT)) +
  geom_point(color = "blue", size = 3, alpha = 0.7) +  
  geom_smooth(method = "lm", color = "red", size = 1.2, linetype = "solid") +  
  labs(title = "Reaction Time by Temporal Demand",
       x = "Temporal Demand",
       y = "Average Block Reaction Time (ms)") +
  theme_minimal() +
  theme(legend.position = "none") 

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 17 rows containing non-finite outside the scale range
## (`stat_smooth()`).

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

# Plotting Residuals and Assumption of Normality

residuals <- resid(LMER_RT_TD)
fitted <- fitted(LMER_RT_TD)

par(mfrow = c(1, 2))
hist(residuals, main = "Histogram of Residuals", xlab = "Residuals", breaks = 20)
qqnorm(residuals)
qqline(residuals, col = "red")

## RT by Performance

LMER_RT_P <- lmer(Avg_Block_RT ~ Performance  + (1 | ID) , data = Merged_Long)

summary(LMER_RT_P)

## Linear mixed model fit by REML ['lmerMod']
## Formula: Avg_Block_RT ~ Performance + (1 | ID)
##    Data: Merged_Long
## 
## REML criterion at convergence: 840.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.5438 -0.4717 -0.1379  0.5252  2.4566 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 46634    216.0   
##  Residual              7691     87.7   
## Number of obs: 68, groups:  ID, 17
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  519.571     70.189   7.402
## Performance   -1.514      3.723  -0.407
## 
## Correlation of Fixed Effects:
##             (Intr)
## Performance -0.648

Plot

ggplot(Merged_Long, aes(x = Performance, y = Avg_Block_RT)) +
  geom_point(color = "blue", size = 3, alpha = 0.7) +  
  geom_smooth(method = "lm", color = "red", size = 1.2, linetype = "solid") +  
  labs(title = "Reaction Time by Performance",
       x = "Performance",
       y = "Average Block Reaction Time (ms)") +
  theme_minimal() +
  theme(legend.position = "none") 

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 17 rows containing non-finite outside the scale range
## (`stat_smooth()`).

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

# Plotting Residuals and Assumption of Normality

residuals <- resid(LMER_RT_P)
fitted <- fitted(LMER_RT_P)

par(mfrow = c(1, 2))
hist(residuals, main = "Histogram of Residuals", xlab = "Residuals", breaks = 20)
qqnorm(residuals)
qqline(residuals, col = "red")

## RT by Effort

LMER_RT_E <- lmer(Avg_Block_RT ~ Effort  + (1 | ID) , data = Merged_Long)

summary(LMER_RT_E)

## Linear mixed model fit by REML ['lmerMod']
## Formula: Avg_Block_RT ~ Effort + (1 | ID)
##    Data: Merged_Long
## 
## REML criterion at convergence: 839.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.5086 -0.5378 -0.1021  0.5137  2.5519 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 46298    215.2   
##  Residual              7638     87.4   
## Number of obs: 68, groups:  ID, 17
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  476.485     61.646   7.729
## Effort         2.637      3.331   0.792
## 
## Correlation of Fixed Effects:
##        (Intr)
## Effort -0.504

Plot

ggplot(Merged_Long, aes(x = Effort, y = Avg_Block_RT)) +
  geom_point(color = "blue", size = 3, alpha = 0.7) +  
  geom_smooth(method = "lm", color = "red", size = 1.2, linetype = "solid") +  
  labs(title = "Reaction Time by Effort",
       x = "Effort",
       y = "Average Block Reaction Time (ms)") +
  theme_minimal() +
  theme(legend.position = "none") 

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 17 rows containing non-finite outside the scale range
## (`stat_smooth()`).

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

# Plotting Residuals and Assumption of Normality

residuals <- resid(LMER_RT_E)
fitted <- fitted(LMER_RT_E)

par(mfrow = c(1, 2))
hist(residuals, main = "Histogram of Residuals", xlab = "Residuals", breaks = 20)
qqnorm(residuals)
qqline(residuals, col = "red")

## RT by Frustration

LMER_RT_F <- lmer(Avg_Block_RT ~ Frustration  + (1 | ID) , data = Merged_Long)

summary(LMER_RT_F)

## Linear mixed model fit by REML ['lmerMod']
## Formula: Avg_Block_RT ~ Frustration + (1 | ID)
##    Data: Merged_Long
## 
## REML criterion at convergence: 840.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.5754 -0.4615 -0.1119  0.5197  2.3887 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 46388    215.38  
##  Residual              7702     87.76  
## Number of obs: 68, groups:  ID, 17
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  490.932     58.678   8.367
## Frustration    1.344      3.250   0.413
## 
## Correlation of Fixed Effects:
##             (Intr)
## Frustration -0.418

Plot

ggplot(Merged_Long, aes(x = Frustration, y = Avg_Block_RT)) +
  geom_point(color = "blue", size = 3, alpha = 0.7) +  
  geom_smooth(method = "lm", color = "red", size = 1.2, linetype = "solid") +  
  labs(title = "Reaction Time by Frustration",
       x = "Frustration",
       y = "Average Block Reaction Time (ms)") +
  theme_minimal() +
  theme(legend.position = "none") 

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 17 rows containing non-finite outside the scale range
## (`stat_smooth()`).

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

# Plotting Residuals and Assumption of Normality

residuals <- resid(LMER_RT_F)
fitted <- fitted(LMER_RT_F)

par(mfrow = c(1, 2))
hist(residuals, main = "Histogram of Residuals", xlab = "Residuals", breaks = 20)
qqnorm(residuals)
qqline(residuals, col = "red")

## RT by Pleasure

LMER_RT_p <- lmer(Avg_Block_RT ~ Pleasure  + (1 | ID) , data = Merged_Long)

summary(LMER_RT_p)

## Linear mixed model fit by REML ['lmerMod']
## Formula: Avg_Block_RT ~ Pleasure + (1 | ID)
##    Data: Merged_Long
## 
## REML criterion at convergence: 826.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.3242 -0.5318 -0.1920  0.4865  2.4499 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 47501    217.95  
##  Residual              7641     87.41  
## Number of obs: 67, groups:  ID, 17
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  440.885     77.088   5.719
## Pleasure       9.848      8.984   1.096
## 
## Correlation of Fixed Effects:
##          (Intr)
## Pleasure -0.715

Plot

ggplot(Merged_Long, aes(x = Pleasure, y = Avg_Block_RT)) +
  geom_point(color = "blue", size = 3, alpha = 0.7) +  
  geom_smooth(method = "lm", color = "red", size = 1.2, linetype = "solid") +  
  labs(title = "Reaction Time by Pleasure",
       x = "Pleasure",
       y = "Average Block Reaction Time (ms)") +
  theme_minimal() +
  theme(legend.position = "none") 

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 18 rows containing non-finite outside the scale range
## (`stat_smooth()`).

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

# Plotting Residuals and Assumption of Normality

residuals <- resid(LMER_RT_p)
fitted <- fitted(LMER_RT_p)

par(mfrow = c(1, 2))
hist(residuals, main = "Histogram of Residuals", xlab = "Residuals", breaks = 20)
qqnorm(residuals)
qqline(residuals, col = "red")

## RT by Arousal

LMER_RT_a <- lmer(Avg_Block_RT ~ Arousal  + (1 | ID)  + (1 | Block), data = Merged_Long)

summary(LMER_RT_a)

## Linear mixed model fit by REML ['lmerMod']
## Formula: Avg_Block_RT ~ Arousal + (1 | ID) + (1 | Block)
##    Data: Merged_Long
## 
## REML criterion at convergence: 826
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.3455 -0.4341 -0.1035  0.3739  2.5101 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 47326.9  217.55  
##  Block    (Intercept)   371.8   19.28  
##  Residual              7394.8   85.99  
## Number of obs: 67, groups:  ID, 17; Block, 4
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  462.976     80.608   5.744
## Arousal        6.698     10.452   0.641
## 
## Correlation of Fixed Effects:
##         (Intr)
## Arousal -0.735

Plot

ggplot(Merged_Long, aes(x = Arousal, y = Avg_Block_RT)) +
  geom_point(color = "blue", size = 3, alpha = 0.7) +  
  geom_smooth(method = "lm", color = "red", size = 1.2, linetype = "solid") +  
  labs(title = "Reaction Time by Arousal",
       x = "Arousal",
       y = "Average Block Reaction Time (ms)") +
  theme_minimal() +
  theme(legend.position = "none")

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 18 rows containing non-finite outside the scale range
## (`stat_smooth()`).

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

# Plotting Residuals and Assumption of Normality

residuals <- resid(LMER_RT_a)
fitted <- fitted(LMER_RT_a)

par(mfrow = c(1, 2))
hist(residuals, main = "Histogram of Residuals", xlab = "Residuals", breaks = 20)
qqnorm(residuals)
qqline(residuals, col = "red")

## RT by BMI

LMER_RT_BMI <- lmer(Response_Time ~ BMI + (1 | ID) + (1 | Block), data = Merged)

Anova(LMER_RT_BMI)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Response_Time
##      Chisq Df Pr(>Chisq)
## BMI 1.0355  1     0.3089

Plot

# Create new dataset with average Block RT variable
Merged_avg <- Merged %>%
  group_by(ID, Block) %>%
  summarise(Avg_Response_Time = mean(Response_Time), .groups = 'drop')

# Join the BMI data to the averaged data
Merged_avg <- Merged_avg %>%
  left_join(Merged %>% select(ID, BMI) %>% distinct(), by = "ID")

# Create the plot using the averaged data
ggplot(Merged_avg, aes(x = BMI, y = Avg_Response_Time)) +
  geom_point(color = "blue", size = 3, alpha = 0.7) +  
  geom_smooth(method = "lm", color = "red", size = 1.2, linetype = "solid") +  
  labs(title = "Reaction Time by BMI",
       x = "BMI",
       y = "Average Block Reaction Time (ms)") +
  theme_minimal() +
  theme(legend.position = "none")

## `geom_smooth()` using formula = 'y ~ x'

# Plotting Residuals and Assumption of Normality

residuals <- resid(LMER_RT_BMI)
fitted <- fitted(LMER_RT_BMI)

par(mfrow = c(1, 2))
hist(residuals, main = "Histogram of Residuals", xlab = "Residuals", breaks = 20)
qqnorm(residuals)
qqline(residuals, col = "red")

RT by Age

LMER_RT_Age <- lmer(Response_Time ~ Age + (1 | ID) + (1 | Block), data = Merged)

Anova(LMER_RT_Age)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Response_Time
##      Chisq Df Pr(>Chisq)
## Age 0.5046  1     0.4775

Plot

# Join the Age data to the averaged data
Merged_avg <- Merged_avg %>%
  left_join(Merged %>% select(ID, Age) %>% distinct(), by = "ID")

# Create the plot using the averaged data
ggplot(Merged_avg, aes(x = Age, y = Avg_Response_Time)) +
  geom_point(color = "blue", size = 3, alpha = 0.7) +  
  geom_smooth(method = "lm", color = "red", size = 1.2, linetype = "solid") +  
  labs(title = "Reaction Time by Age",
       x = "Age",
       y = "Average Block Reaction Time (ms)") +
  theme_minimal() +
  theme(legend.position = "none")

## `geom_smooth()` using formula = 'y ~ x'

# Plotting Residuals and Assumption of Normality

residuals <- resid(LMER_RT_Age)
fitted <- fitted(LMER_RT_Age)

par(mfrow = c(1, 2))
hist(residuals, main = "Histogram of Residuals", xlab = "Residuals", breaks = 20)
qqnorm(residuals)
qqline(residuals, col = "red")

RT by Gender

LMER_RT_Gender <- lmer(Response_Time ~ Gender + (1 | ID) + (1 | Block), data = Merged)

Anova(LMER_RT_Gender)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Response_Time
##         Chisq Df Pr(>Chisq)  
## Gender 3.5326  1    0.06017 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Plot

# Join the Gender data to the averaged data
Merged_avg <- Merged_avg %>%
  left_join(Merged %>% select(ID, Gender) %>% distinct(), by = "ID")

# Ensure Gender is a factor variable
Merged_avg$Gender <- as.factor(Merged_avg$Gender)

# Create the plot using the averaged data as a boxplot
ggplot(Merged_avg, aes(x = Gender, y = Avg_Response_Time)) +
  geom_boxplot(fill = "skyblue", color = "black", alpha = 0.7) +  
  labs(title = "Reaction Time by Gender",
       x = "Gender",
       y = "Average Block Reaction Time (ms)") +
  theme_minimal() +
  theme(legend.position = "none")

RT by Handedness

LMER_RT_Hand <- lmer(Response_Time ~ Footedness + (1 | ID) + (1 | Block), data = Merged)

Anova(LMER_RT_Hand)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Response_Time
##            Chisq Df Pr(>Chisq)
## Footedness 0.039  1     0.8435

Plot

# Join the Hand data to the averaged data
Merged_avg <- Merged_avg %>%
  left_join(Merged %>% select(ID, Footedness) %>% distinct(), by = "ID")

# Ensure Hand is a factor variable
Merged_avg$Gender <- as.factor(Merged_avg$Gender)

# Create the plot using the averaged data as a boxplot
ggplot(Merged_avg, aes(x = Gender, y = Avg_Response_Time)) +
  geom_boxplot(fill = "skyblue", color = "black", alpha = 0.7) +  
  labs(title = "Reaction Time by Handedness",
       x = "Handedness",
       y = "Average Block Reaction Time (ms)") +
  theme_minimal() +
  theme(legend.position = "none")

RT by Smoking

LMER_RT_Smoke <- lmer(Response_Time ~ Smoking + (1 | ID) + (1 | Block), data = Merged)

Anova(LMER_RT_Smoke)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Response_Time
##          Chisq Df Pr(>Chisq)  
## Smoking 4.2566  1     0.0391 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Plot

# Join the Smoking data to the averaged data
Merged_avg <- Merged_avg %>%
  left_join(Merged %>% select(ID, Smoking) %>% distinct(), by = "ID")

# Ensure Smoking is a factor variable
Merged_avg$Smoking <- as.factor(Merged_avg$Smoking)

# Create the plot using the averaged data as a boxplot
ggplot(Merged_avg, aes(x = Smoking, y = Avg_Response_Time)) +
  geom_boxplot(fill = "skyblue", color = "black", alpha = 0.7) +  
  labs(title = "Reaction Time by Smoking",
       x = "Smoking",
       y = "Average Block Reaction Time (ms)") +
  theme_minimal() +
  theme(legend.position = "none")

RT by Alcohol

LMER_RT_Alc <- lmer(Response_Time ~ Alcohol + (1 | ID) + (1 | Block), data = Merged)

Anova(LMER_RT_Alc)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Response_Time
##         Chisq Df Pr(>Chisq)
## Alcohol  0.14  1     0.7083

Plot

# Join the Alcohol data to the averaged data
Merged_avg <- Merged_avg %>%
  left_join(Merged %>% select(ID, Alcohol) %>% distinct(), by = "ID")

# Ensure Alcohol is a factor variable
Merged_avg$Alcohol <- as.factor(Merged_avg$Alcohol)

# Create the plot using the averaged data as a boxplot
ggplot(Merged_avg, aes(x = Alcohol, y = Avg_Response_Time)) +
  geom_boxplot(fill = "skyblue", color = "black", alpha = 0.7) +  
  labs(title = "Reaction Time by Alcohol Consumption (last 24h)",
       x = "Alcohol Consumption",
       y = "Average Block Reaction Time (ms)") +
  theme_minimal() +
  theme(legend.position = "none")