Multimodal DS_DSP Analysis: Behavioural, DTW, ERD/S, Models and Post-Hoc
Victoria Lakomski, Johannes Pannermayr
2024-05-15
##Load relevant Libraries
library(lme4)
library(readxl)
library(haven)
library(effects)
library(lattice)
library(car)
library(ggplot2)
library(knitr)
library(reshape2)
library(dplyr)
library(forcats)
library(DHARMa)
library(Hmisc)
library(phia)
library(lsmeans)
library(emmeans)
library(multcomp)
library(plotly)
library(lmerTest)
library(tinytex)
library(tidyverse)
library(devtools)
library(brms)
library(janitor)
library(magrittr)
library(mascutils)
library(ggthemes)
library(openxlsx)##Load dataset, remove irrelevant variables and replace participant 1 with 13
df <- read_excel("~/uni/data thesis/df_P23678.xlsx")
# remove rows not used for analysis (mainly cue presentation and sequence markers)
df <- df[!(df$procedure %in% c("cueprocedure", "fixatieprocedure", "toetsenbordprocedure", "cuelist2", "fixprocedure", "seq1", "seq2", "sequentieprocedure", "feedbackprocedure", "nogoprocedure")), ]
df <- subset(df, select = -c(cue.OnsetTime, cue.OnsetDelay))
#remove participant 1
df<- df %>%
filter(subject != 1)
#rename participant 13 into 1
df <- df %>%
mutate(subject = ifelse(subject == 13, 1, subject))##arrange so that sessions are in correct order
df <- df %>%
arrange(subject, session)#add count variable to include the trial number
df <- df %>%
group_by(subject) %>%
mutate(trial_nr = row_number() %% 6 == 0) %>%
mutate(trial_nr = cumsum(trial_nr)) %>% ungroup()##remove incorrect As soon as 1 step is incorrect within a sequence we reject the whole thing
num_blocks <- nrow(df)/ 6
rows_to_remove <- c()
for(i in 0:(num_blocks - 1)) {
block_indices <- (1:6) + (i * 6)
if(any(df$feedback.ACC[block_indices] == 0)){
rows_to_remove <- c(rows_to_remove, block_indices)}}
df <- df[-rows_to_remove, ]
df <- df %>%
select(-procedure, -feedback.ACC)
df$feedback.RT <- as.numeric(df$feedback.RT)##now we add a count variable to count the correct sequence and rename the sub trial variable when you look at the dataset now, the correct_trial and tiral_nr variable seem to not overlap because the trial_nr variable starts with 0, but they do always overlap on the 6th step which is the crucial one because we create the averages there. When we create the averages dataframe this issue does not matter anymore.
df <- df %>%
mutate(correct_trial = ((row_number() - 1) %/% 6) + 1)
df <- df %>%
rename(stepNumber = `sub trial number`)
df = df %>% select(subject, session, stepNumber, feedback.RT, h, trial_nr, correct_trial)##calculate the ave RT and sum RT we need to make a variable called sequence_id to give each sequence a unique number. Now we have repeating numbers because our trial counts reset per participant.
df <- df %>%
mutate(sequence_id = rep(1:(n() %/% 6), each = 6)) %>%
group_by(sequence_id) %>%
mutate(aveRT= mean(feedback.RT), sumRT = sum(feedback.RT)) %>% ungroup()##removing outliers using IQR Do this when the EEG data is bound. 258 trials (43 sequences) removed
Q1 <- quantile(df$aveRT, 0.25)
Q3 <- quantile(df$aveRT, 0.75)
IQR <- Q3 - Q1
upper_bound <- Q3 + 1.5 * IQR
lower_bound <- Q1 - 1.5 * IQR
df <- df %>%
filter(aveRT < upper_bound)##make a seperate dataset which does not include the steps for learning curves
averages = df %>% filter(stepNumber == 6)##save datasets
setwd("~/uni/data thesis")
FINAL_DATA <- FINAL_DATA %>%
rename(step = `stepNumber`)
FINAL_DATA = FINAL_DATA %>% rename(block = `session`)
write.xlsx(averages, "averages.xlsx")
write.xlsx(FINAL_DATA, "FINAL_DATA.xlsx")#This is the EEG data DTW step This chunk needs to be run for each participant and block. You can change the participant number and block number at the top of the chunk to make sure the resulting file is named correctly. The code will go through all files (= trials) in the folder and store them in a list. Then it will go through this list and apply the DTW mechanism and later the ERD/S calculation. The result is one file per participant and block containing all trials on step level basis. The chunk is currently deactivated.
setwd("C:/Users/Alex/Desktop/eeg data/participant5/theta/block4")
#set your wd to folder containing block trials
##read all excel filed from folder
excel_files <- list.files(pattern = "\\.xlsx$")
data_frames <- list()
Participant <- 5
Block <- 4
file_name <- paste0("P", Participant, "_", Block, ".xlsx")
# Loop through each Excel file and import it into R
for (file in excel_files) {
# Read the Excel file into a data frame
data <- read_excel(file)
# Store the data frame in the list
data_frames[[file]] <- data
}
##select only channels c3. c4 and cz
select_columns <- function(df) {
df <- df[, c('C3', 'C4', 'Cz')]
return(df)
}
# Apply the function to each data frame in the list
data_frames <- lapply(data_frames, select_columns)
##remove the first 101 rows of all datasets because that is before the 27 marker (200ms) and save in a seperate list
remove_first_rows <- function(df) {
first_rows <- df[1:101, ]
df <- df[-(1:101), ]
return(list(first_rows = first_rows, modified_df = df))
}
# Apply the function to each data frame in the list
baseline <- lapply(data_frames, function(df) remove_first_rows(df)$first_rows)
modified_data <- lapply(data_frames, function(df) remove_first_rows(df)$modified_df)
#remove last 100ms of all datasets which is 51 rows
remove_last_rows <- function(df) {
df <- df[1:(nrow(df) - 51), ]
return(df)
}
# Apply the function to each data frame in the modified_data list
final_modified_data <- lapply(modified_data, remove_last_rows)
# Define a function to downsize the dataset
downsize_dataset <- function(df) {
# Calculate the number of rows in the dataset
num_rows <- nrow(df)
# Define the number of bins (always use 300 bins)
num_bins <- 300
# Calculate the number of rows per bin
rows_per_bin <- floor(num_rows / num_bins)
# Calculate the number of excess rows
excess_rows <- num_rows %% num_bins
# Calculate the bin width
bin_width <- rep(rows_per_bin, num_bins)
# Distribute excess rows evenly among bins
if (excess_rows > 0) {
bin_width[1:excess_rows] <- bin_width[1:excess_rows] + 1
}
# Create bin indices
bin_indices <- rep(1:num_bins, times = bin_width)
# Calculate the mean of the values within each bin for each column separately
df_binned <- df %>%
group_by(bin = bin_indices) %>%
summarise(mean_C3 = mean(C3),
mean_C4 = mean(C4),
mean_Cz = mean(Cz), .groups = "drop")
return(df_binned)
}
# Apply the function to each data frame in the final_modified_data list
downsized_data <- lapply(final_modified_data, downsize_dataset)
##get mean for baseline in each dataset
baseline_mean = function(df)
colMeans(df, na.rm = TRUE)
mean_baseline <- lapply(baseline, baseline_mean)
##make step dataset for each trial
calculate_averages <- function(df) {
# Add a new variable called "step" indicating the group number (each group consists of 50 bins)
df <- df %>%
mutate(step = ceiling(bin / 50))
# Calculate the average for each group of 50 bins
averaged_df <- df %>%
group_by(step) %>%
summarise(C3 = mean(mean_C3),
C4 = mean(mean_C4),
Cz = mean(mean_Cz), .groups = "drop")
return(averaged_df)
}
# Apply the function to each dataframe in final_modified_data list
averaged_trials <- lapply(downsized_data, calculate_averages)
##add baseline value to each step dataset (per trial)
# Define a function to add the Step column with value 0 to each dataset in baseline_list
add_step_column <- function(df) {
# Add the Step column with value 0
df$step <- 0
return(df)
}
# Apply the function to each dataset in the baseline_list
mean_baseline <- lapply(mean_baseline, add_step_column)
for (i in seq_along(mean_baseline)) {
# Add the observation to the corresponding dataset in list2
averaged_trials[[i]] <- rbind(mean_baseline[[i]], averaged_trials[[i]])
}
##change step column to chracter so the baseline calculation is not iterated over it
##convert time to numeric
convert_step_to_character <- function(df) {
df$step <- as.character(df$step)
return(df)
}
# Apply the function to each data frame in the list
averaged_trials <- lapply(averaged_trials, convert_step_to_character)
save_file = averaged_trials
##apply ERD/S formula to the datasets
# Iterate over each dataset in the list
for (i in seq_along(averaged_trials)) {
# Get the current dataset
df <- averaged_trials[[i]]
# Get the baseline values from row 1
baseline <- df[1, ]
# Iterate over each column
for (col_name in colnames(df)) {
# Skip the column if it's not numeric or if it's the first column (which is the baseline)
if (!is.numeric(df[[col_name]])) {
next
}
# Apply the formula to calculate ERDS
df[[col_name]] <- (df[[col_name]] - baseline[[col_name]]) / baseline[[col_name]] * 100
}
# Assign the updated dataset back to the list
averaged_trials[[i]] <- df
}
##change step to numeric again
convert_step_to_numeric <- function(df) {
df$step <- as.numeric(df$step)
return(df)
}
# Apply the function to each data frame in the list
averaged_trials <- lapply(averaged_trials, convert_step_to_numeric)
##add participant number, block number, trial number
##participant and block number needs to be changed every time
for (i in seq_along(averaged_trials)) {
# Add a new column "subject" with the value 5 to the current dataset
averaged_trials[[i]]$subject <- Participant
}
for (i in seq_along(averaged_trials)) {
# Add a new column "subject" with the value 5 to the current dataset
averaged_trials[[i]]$block <- Block
}
for (i in seq_along(averaged_trials)) {
# Add a new column "trial" with the corresponding trial number to the current dataset
averaged_trials[[i]]$trial <- i
}
#combine all datasets and remove baseline
df2 <- do.call(rbind, averaged_trials)
df2 = df2 %>% filter(step != 0)
df2 = df2 %>%
rownames_to_column() %>%
select(-rowname)
# save this file
library(openxlsx)
write.xlsx(df2, file_name)##Binding all files The next code goes through all the other files previously made and binds them into one large file, just as the behavioural data is. Also currently deactivated. Make sure all the previous files are in one folder and set WD.
#RT Models and plots over all sessions/blocks
setwd("C:/Users/Alex/Desktop/eeg data")
averages <- read_excel("C:/Users/Alex/Desktop/eeg data/averages.xlsx")
#make factors
averages$subject = factor(averages$subject)
averages$block = factor(averages$block, ordered = TRUE, levels=c('1', '2', '3', '4', '5', '6'))
#RT model
m.RT = lmer(aveRT ~ block + (1|subject), data = averages, REML = FALSE)
#plot residuals
residuals <- resid(m.RT)
fitted <- fitted(m.RT)
par(mfrow = c(1, 2))
hist(residuals, main = "Histogram of Residuals", xlab = "Residuals", breaks = 20)
qqnorm(residuals)
qqline(residuals, col = "red")
Anova(m.RT)summary(m.RT, ddf = "Satterthwaite")## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: aveRT ~ block + (1 | subject)
## Data: averages
##
## AIC BIC logLik deviance df.resid
## 34133.9 34181.4 -17058.9 34117.9 2816
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.5527 -0.5805 -0.0780 0.4504 5.2548
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 18777 137.0
## Residual 10074 100.4
## Number of obs: 2824, groups: subject, 12
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 360.200 39.603 11.998 9.095 9.89e-07 ***
## block.L -152.077 4.756 2812.075 -31.975 < 2e-16 ***
## block.Q 137.232 4.726 2812.030 29.037 < 2e-16 ***
## block.C 5.067 4.617 2812.020 1.097 0.27261
## block^4 31.705 4.554 2812.006 6.961 4.17e-12 ***
## block^5 14.616 4.554 2812.001 3.209 0.00135 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) blck.L blck.Q blck.C blck^4
## block.L -0.001
## block.Q 0.003 -0.048
## block.C -0.001 0.058 -0.026
## block^4 0.001 -0.018 0.022 -0.008
## block^5 0.000 0.013 -0.008 0.000 0.004#posthoc
emm.RT = emmeans(m.RT, pairwise ~ block)
print(emm.RT)## $emmeans
## block emmean SE df lower.CL upper.CL
## 1 529 41.6 13.4 440 619
## 2 389 41.6 13.4 299 479
## 3 323 41.6 13.4 233 412
## 4 302 41.6 13.4 212 391
## 5 265 41.6 13.4 176 355
## 6 353 41.6 13.4 263 443
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## block1 - block2 140.1 6.73 2817 20.818 <.0001
## block1 - block3 206.4 6.76 2817 30.518 <.0001
## block1 - block4 227.3 6.79 2817 33.492 <.0001
## block1 - block5 263.7 6.76 2817 39.007 <.0001
## block1 - block6 176.1 6.85 2817 25.710 <.0001
## block2 - block3 66.3 6.40 2817 10.347 <.0001
## block2 - block4 87.2 6.43 2817 13.562 <.0001
## block2 - block5 123.6 6.40 2817 19.303 <.0001
## block2 - block6 36.0 6.50 2817 5.548 <.0001
## block3 - block4 21.0 6.46 2817 3.246 0.0150
## block3 - block5 57.3 6.43 2817 8.915 <.0001
## block3 - block6 -30.2 6.52 2817 -4.632 0.0001
## block4 - block5 36.3 6.45 2817 5.634 <.0001
## block4 - block6 -51.2 6.55 2817 -7.818 <.0001
## block5 - block6 -87.5 6.52 2817 -13.426 <.0001
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 6 estimates#Effects
ae.m.RT <- allEffects(m.RT)
df.ae.m.RT <- as.data.frame(ae.m.RT[[1]])
#Proper CIs
#change conf interval to 83%, match p = .05
df.ae.m.RT$l83 <- df.ae.m.RT$fit - 1.3722 * df.ae.m.RT$se
df.ae.m.RT$u83 <- df.ae.m.RT$fit + 1.3722 * df.ae.m.RT$se
# Effects Plot
plot.m.RT <- ggplot(df.ae.m.RT, 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() +
#geom_path(aes(x=session, y=fit, group=session, colour=session)) +
#geom_ribbon(aes(ymin=l83, ymax=u83, group=session, fill=session), alpha=0.2) +
ylab("RT (ms)")+
xlab("Block")+
ggtitle("Learning (1-4) and Test (5-6) Phase RT")+
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.m.RT)## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
#Concatenation
setwd("C:/Users/Alex/Desktop/eeg data")
FINAL_DATA <- read_excel("C:/Users/Alex/Desktop/eeg data/FINAL_DATA.xlsx")
#make factors
FINAL_DATA$subject = factor(FINAL_DATA$subject)
FINAL_DATA$block = factor(FINAL_DATA$block, ordered = TRUE, levels=c('1', '2', '3', '4', '5', '6'))
FINAL_DATA$step= factor(FINAL_DATA$step, ordered = TRUE, levels=c('1' ,'2', '3', '4', '5', '6'))
#RT model
m.conRT = lmer(feedback.RT ~ step * block + (1|subject), data = FINAL_DATA, REML = FALSE)
#plot residuals
residuals2 <- resid(m.conRT)
fitted2 <- fitted(m.conRT)
par(mfrow = c(1, 2))
hist(residuals2, main = "Histogram of Residuals", xlab = "Residuals", breaks = 20)
qqnorm(residuals2)
qqline(residuals2, col = "red")
anova(m.conRT)summary(m.conRT, dff = "Satterthwaite")## Warning in summary.merMod(as(object, "lmerMod"), ...): additional arguments
## ignored## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: feedback.RT ~ step * block + (1 | subject)
## Data: FINAL_DATA
##
## AIC BIC logLik deviance df.resid
## 228859.0 229153.0 -114391.5 228783.0 16906
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.2269 -0.5008 -0.1123 0.3313 16.3526
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 14744 121.4
## Residual 42618 206.4
## Number of obs: 16944, groups: subject, 12
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 357.8493 35.0890 12.0002 10.198 2.89e-07 ***
## step.L -239.8845 3.8908 16932.0001 -61.655 < 2e-16 ***
## step.Q 237.4983 3.8908 16932.0001 61.042 < 2e-16 ***
## step.C -146.1383 3.8908 16932.0001 -37.560 < 2e-16 ***
## step^4 63.1122 3.8908 16932.0001 16.221 < 2e-16 ***
## step^5 -35.3115 3.8908 16932.0001 -9.076 < 2e-16 ***
## block.L -82.0607 3.9751 16932.3105 -20.644 < 2e-16 ***
## block.Q 7.6317 3.9579 16932.1815 1.928 0.053848 .
## block.C 45.3485 3.8747 16932.1010 11.704 < 2e-16 ***
## block^4 -9.3205 3.8291 16932.0454 -2.434 0.014939 *
## block^5 8.5690 3.8337 16932.0349 2.235 0.025418 *
## step.L:block.L 42.3509 9.7158 16932.0001 4.359 1.31e-05 ***
## step.Q:block.L -13.5569 9.7158 16932.0001 -1.395 0.162932
## step.C:block.L 9.9448 9.7158 16932.0001 1.024 0.306054
## step^4:block.L 1.6275 9.7158 16932.0001 0.168 0.866967
## step^5:block.L 12.9141 9.7158 16932.0001 1.329 0.183806
## step.L:block.Q 8.2767 9.6820 16932.0001 0.855 0.392648
## step.Q:block.Q -16.1845 9.6820 16932.0001 -1.672 0.094622 .
## step.C:block.Q -7.5802 9.6820 16932.0001 -0.783 0.433691
## step^4:block.Q -4.7546 9.6820 16932.0001 -0.491 0.623379
## step^5:block.Q 22.8010 9.6820 16932.0001 2.355 0.018535 *
## step.L:block.C -36.3369 9.4842 16932.0001 -3.831 0.000128 ***
## step.Q:block.C 0.8732 9.4842 16932.0001 0.092 0.926641
## step.C:block.C -17.7732 9.4842 16932.0001 -1.874 0.060951 .
## step^4:block.C -1.2547 9.4842 16932.0001 -0.132 0.894750
## step^5:block.C -11.0398 9.4842 16932.0001 -1.164 0.244431
## step.L:block^4 -9.7621 9.3763 16932.0001 -1.041 0.297820
## step.Q:block^4 11.1006 9.3763 16932.0001 1.184 0.236467
## step.C:block^4 -4.4267 9.3763 16932.0001 -0.472 0.636845
## step^4:block^4 -11.1748 9.3763 16932.0001 -1.192 0.233348
## step^5:block^4 -9.5678 9.3763 16932.0001 -1.020 0.307539
## step.L:block^5 9.5879 9.3881 16932.0001 1.021 0.307138
## step.Q:block^5 -7.4790 9.3881 16932.0001 -0.797 0.425666
## step.C:block^5 15.1025 9.3881 16932.0001 1.609 0.107705
## step^4:block^5 -6.0371 9.3881 16932.0001 -0.643 0.520192
## step^5:block^5 5.1378 9.3881 16932.0001 0.547 0.584202
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Correlation matrix not shown by default, as p = 36 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need itemm.conRT = emmeans(m.conRT, pairwise ~ step|block)## Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
## To enable adjustments, add the argument 'pbkrtest.limit = 16944' (or larger)
## [or, globally, 'set emm_options(pbkrtest.limit = 16944)' 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 = 16944' (or larger)
## [or, globally, 'set emm_options(lmerTest.limit = 16944)' or larger];
## but be warned that this may result in large computation time and memory use.print(emm.conRT)## $emmeans
## block = 1:
## step emmean SE df asymp.LCL asymp.UCL
## 1 741 36.5 Inf 669 812
## 2 334 36.5 Inf 263 406
## 3 312 36.5 Inf 241 384
## 4 307 36.5 Inf 235 378
## 5 331 36.5 Inf 260 403
## 6 326 36.5 Inf 254 397
##
## block = 2:
## step emmean SE df asymp.LCL asymp.UCL
## 1 777 36.3 Inf 706 848
## 2 358 36.3 Inf 287 429
## 3 353 36.3 Inf 282 424
## 4 328 36.3 Inf 257 399
## 5 351 36.3 Inf 280 422
## 6 341 36.3 Inf 270 412
##
## block = 3:
## step emmean SE df asymp.LCL asymp.UCL
## 1 745 36.3 Inf 673 816
## 2 299 36.3 Inf 228 371
## 3 302 36.3 Inf 231 373
## 4 262 36.3 Inf 190 333
## 5 310 36.3 Inf 239 381
## 6 296 36.3 Inf 225 367
##
## block = 4:
## step emmean SE df asymp.LCL asymp.UCL
## 1 661 36.3 Inf 590 732
## 2 275 36.3 Inf 204 346
## 3 264 36.3 Inf 192 335
## 4 235 36.3 Inf 164 306
## 5 273 36.3 Inf 202 344
## 6 291 36.3 Inf 220 362
##
## block = 5:
## step emmean SE df asymp.LCL asymp.UCL
## 1 619 36.3 Inf 548 691
## 2 235 36.3 Inf 164 307
## 3 233 36.3 Inf 162 304
## 4 235 36.3 Inf 164 306
## 5 245 36.3 Inf 173 316
## 6 272 36.3 Inf 201 343
##
## block = 6:
## step emmean SE df asymp.LCL asymp.UCL
## 1 653 36.4 Inf 582 724
## 2 267 36.4 Inf 195 338
## 3 249 36.4 Inf 177 320
## 4 259 36.4 Inf 187 330
## 5 276 36.4 Inf 205 347
## 6 269 36.4 Inf 198 340
##
## Degrees-of-freedom method: asymptotic
## Confidence level used: 0.95
##
## $contrasts
## block = 1:
## contrast estimate SE df z.ratio p.value
## step1 - step2 406.40 14.2 Inf 28.596 <.0001
## step1 - step3 428.35 14.2 Inf 30.140 <.0001
## step1 - step4 433.79 14.2 Inf 30.522 <.0001
## step1 - step5 409.40 14.2 Inf 28.807 <.0001
## step1 - step6 414.75 14.2 Inf 29.183 <.0001
## step2 - step3 21.95 14.2 Inf 1.544 0.6355
## step2 - step4 27.38 14.2 Inf 1.927 0.3856
## step2 - step5 3.00 14.2 Inf 0.211 0.9999
## step2 - step6 8.35 14.2 Inf 0.587 0.9919
## step3 - step4 5.44 14.2 Inf 0.383 0.9989
## step3 - step5 -18.95 14.2 Inf -1.333 0.7666
## step3 - step6 -13.60 14.2 Inf -0.957 0.9314
## step4 - step5 -24.38 14.2 Inf -1.716 0.5212
## step4 - step6 -19.04 14.2 Inf -1.340 0.7629
## step5 - step6 5.35 14.2 Inf 0.376 0.9990
##
## block = 2:
## contrast estimate SE df z.ratio p.value
## step1 - step2 419.15 13.1 Inf 32.006 <.0001
## step1 - step3 424.29 13.1 Inf 32.399 <.0001
## step1 - step4 449.84 13.1 Inf 34.350 <.0001
## step1 - step5 426.81 13.1 Inf 32.592 <.0001
## step1 - step6 436.03 13.1 Inf 33.295 <.0001
## step2 - step3 5.14 13.1 Inf 0.392 0.9988
## step2 - step4 30.69 13.1 Inf 2.343 0.1769
## step2 - step5 7.66 13.1 Inf 0.585 0.9920
## step2 - step6 16.88 13.1 Inf 1.289 0.7915
## step3 - step4 25.55 13.1 Inf 1.951 0.3711
## step3 - step5 2.52 13.1 Inf 0.193 1.0000
## step3 - step6 11.74 13.1 Inf 0.896 0.9475
## step4 - step5 -23.02 13.1 Inf -1.758 0.4931
## step4 - step6 -13.81 13.1 Inf -1.055 0.8992
## step5 - step6 9.21 13.1 Inf 0.704 0.9816
##
## block = 3:
## contrast estimate SE df z.ratio p.value
## step1 - step2 445.11 13.3 Inf 33.507 <.0001
## step1 - step3 442.66 13.3 Inf 33.322 <.0001
## step1 - step4 483.04 13.3 Inf 36.362 <.0001
## step1 - step5 434.75 13.3 Inf 32.726 <.0001
## step1 - step6 448.47 13.3 Inf 33.759 <.0001
## step2 - step3 -2.45 13.3 Inf -0.184 1.0000
## step2 - step4 37.93 13.3 Inf 2.855 0.0493
## step2 - step5 -10.36 13.3 Inf -0.780 0.9709
## step2 - step6 3.36 13.3 Inf 0.253 0.9999
## step3 - step4 40.38 13.3 Inf 3.039 0.0286
## step3 - step5 -7.92 13.3 Inf -0.596 0.9913
## step3 - step6 5.81 13.3 Inf 0.437 0.9980
## step4 - step5 -48.29 13.3 Inf -3.635 0.0038
## step4 - step6 -34.57 13.3 Inf -2.602 0.0966
## step5 - step6 13.72 13.3 Inf 1.033 0.9069
##
## block = 4:
## contrast estimate SE df z.ratio p.value
## step1 - step2 386.38 13.3 Inf 29.025 <.0001
## step1 - step3 397.65 13.3 Inf 29.872 <.0001
## step1 - step4 426.28 13.3 Inf 32.023 <.0001
## step1 - step5 388.12 13.3 Inf 29.156 <.0001
## step1 - step6 370.15 13.3 Inf 27.806 <.0001
## step2 - step3 11.27 13.3 Inf 0.847 0.9587
## step2 - step4 39.90 13.3 Inf 2.998 0.0325
## step2 - step5 1.74 13.3 Inf 0.131 1.0000
## step2 - step6 -16.23 13.3 Inf -1.219 0.8276
## step3 - step4 28.63 13.3 Inf 2.151 0.2612
## step3 - step5 -9.53 13.3 Inf -0.716 0.9801
## step3 - step6 -27.50 13.3 Inf -2.066 0.3053
## step4 - step5 -38.17 13.3 Inf -2.867 0.0476
## step4 - step6 -56.14 13.3 Inf -4.217 0.0004
## step5 - step6 -17.97 13.3 Inf -1.350 0.7570
##
## block = 5:
## contrast estimate SE df z.ratio p.value
## step1 - step2 384.05 13.2 Inf 29.030 <.0001
## step1 - step3 386.63 13.2 Inf 29.225 <.0001
## step1 - step4 384.83 13.2 Inf 29.089 <.0001
## step1 - step5 374.99 13.2 Inf 28.345 <.0001
## step1 - step6 347.37 13.2 Inf 26.257 <.0001
## step2 - step3 2.58 13.2 Inf 0.195 1.0000
## step2 - step4 0.78 13.2 Inf 0.059 1.0000
## step2 - step5 -9.06 13.2 Inf -0.685 0.9837
## step2 - step6 -36.69 13.2 Inf -2.773 0.0618
## step3 - step4 -1.80 13.2 Inf -0.136 1.0000
## step3 - step5 -11.64 13.2 Inf -0.880 0.9514
## step3 - step6 -39.26 13.2 Inf -2.968 0.0355
## step4 - step5 -9.84 13.2 Inf -0.744 0.9764
## step4 - step6 -37.47 13.2 Inf -2.832 0.0526
## step5 - step6 -27.63 13.2 Inf -2.088 0.2934
##
## block = 6:
## contrast estimate SE df z.ratio p.value
## step1 - step2 386.43 13.7 Inf 28.202 <.0001
## step1 - step3 404.36 13.7 Inf 29.511 <.0001
## step1 - step4 394.31 13.7 Inf 28.778 <.0001
## step1 - step5 377.19 13.7 Inf 27.528 <.0001
## step1 - step6 383.78 13.7 Inf 28.009 <.0001
## step2 - step3 17.94 13.7 Inf 1.309 0.7801
## step2 - step4 7.89 13.7 Inf 0.576 0.9926
## step2 - step5 -9.23 13.7 Inf -0.674 0.9848
## step2 - step6 -2.64 13.7 Inf -0.193 1.0000
## step3 - step4 -10.05 13.7 Inf -0.733 0.9778
## step3 - step5 -27.17 13.7 Inf -1.983 0.3521
## step3 - step6 -20.58 13.7 Inf -1.502 0.6630
## step4 - step5 -17.12 13.7 Inf -1.250 0.8122
## step4 - step6 -10.53 13.7 Inf -0.769 0.9727
## step5 - step6 6.59 13.7 Inf 0.481 0.9968
##
## Degrees-of-freedom method: asymptotic
## P value adjustment: tukey method for comparing a family of 6 estimates#concat effects and plots
ae.m.conRT <- allEffects(m.conRT)
df.ae.m.conRT <- as.data.frame(ae.m.conRT[[1]])
#Proper CIs
#change conf interval to 83%, match p = .05
df.ae.m.conRT$l83 <- df.ae.m.conRT$fit - 1.3722 * df.ae.m.conRT$se
df.ae.m.conRT$u83 <- df.ae.m.conRT$fit + 1.3722 * df.ae.m.conRT$se
# Effects Plot
plot.m.conRT <- ggplot(df.ae.m.conRT, aes(x=step, y=fit, ymin=l83, ymax=u83))+
geom_point(aes(color = block), size=3) +
geom_path(aes(x=step, y=fit, group=block, colour=block)) +
geom_ribbon(aes(ymin=l83, ymax=u83, group=block, fill=block), alpha=0.07) +
ylab("RT (ms)")+
xlab("Step Position")+
ggtitle("Behavioual Data Concatenation")+
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))
print(plot.m.conRT)
plot.m.conRT2 <- ggplot(df.ae.m.conRT, aes(x=step, y=fit, ymin=l83, ymax=u83))+
geom_point(aes(color = block), size=3) +
geom_path(aes(x=step, y=fit, group=block, colour=block)) +
geom_ribbon(aes(ymin=l83, ymax=u83, group=block, fill=block), alpha=0.2) +
ylab("RT (ms)")+
xlab("Step Position")+
ggtitle("Concatenation RT")+
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))+
facet_wrap(~block)
print(plot.m.conRT2)
#Theta Analysis ##Making models with C3, C4 Cz to determine AIC and BIC for futher continuation
FINAL_DATA$subject = factor(FINAL_DATA$subject)
FINAL_DATA$block = factor(FINAL_DATA$block, ordered = TRUE, levels=c('1', '2', '3', '4', '5', '6'))
FINAL_DATA$step = factor(FINAL_DATA$step, ordered = TRUE, levels=c('1', '2', '3', '4', '5', '6'))
m.thetaC3 <- lmer(C3 ~ block + (1|subject), data=FINAL_DATA, REML = FALSE)
m.thetaC4 <- lmer(C4 ~ block + (1|subject), data=FINAL_DATA, REML = FALSE)
m.thetaCz <- lmer(Cz ~ block + (1|subject), data=FINAL_DATA, REML = FALSE)
AIC (m.thetaC3, m.thetaC4, m.thetaCz)BIC (m.thetaC3, m.thetaC4, m.thetaCz)#Continue with Cz bc highest AIC and BIC
anova(m.thetaCz)summary(m.thetaCz)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: Cz ~ block + (1 | subject)
## Data: FINAL_DATA
##
## AIC BIC logLik deviance df.resid
## 302759.2 302821.1 -151371.6 302743.2 16936
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.584 -0.201 -0.105 0.022 36.922
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 449666 670.6
## Residual 3354082 1831.4
## Number of obs: 16944, groups: subject, 12
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 508.24 194.09 12.01 2.619 0.02244 *
## block.L -103.49 35.26 16932.80 -2.935 0.00334 **
## block.Q 56.95 35.11 16932.47 1.622 0.10484
## block.C 241.74 34.37 16932.27 7.033 2.1e-12 ***
## block^4 27.64 33.97 16932.12 0.814 0.41590
## block^5 -88.55 34.01 16932.10 -2.604 0.00923 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) blck.L blck.Q blck.C blck^4
## block.L -0.001
## block.Q 0.003 -0.028
## block.C -0.001 0.051 -0.015
## block^4 0.002 -0.013 0.019 0.000
## block^5 0.000 0.012 -0.003 -0.002 0.004emm.thetaCz <- emmeans(m.thetaCz, pairwise ~ block)## Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
## To enable adjustments, add the argument 'pbkrtest.limit = 16944' (or larger)
## [or, globally, 'set emm_options(pbkrtest.limit = 16944)' 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 = 16944' (or larger)
## [or, globally, 'set emm_options(lmerTest.limit = 16944)' or larger];
## but be warned that this may result in large computation time and memory use.print(emm.thetaCz)## $emmeans
## block emmean SE df asymp.LCL asymp.UCL
## 1 522 197 Inf 135.8 908
## 2 622 196 Inf 236.6 1007
## 3 634 197 Inf 248.8 1019
## 4 354 197 Inf -31.7 739
## 5 351 197 Inf -34.2 736
## 6 567 197 Inf 181.6 953
##
## Degrees-of-freedom method: asymptotic
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df z.ratio p.value
## block1 - block2 -99.84 49.6 Inf -2.015 0.3338
## block1 - block3 -112.19 49.9 Inf -2.248 0.2160
## block1 - block4 168.25 50.0 Inf 3.367 0.0099
## block1 - block5 170.85 49.8 Inf 3.428 0.0080
## block1 - block6 -45.34 50.7 Inf -0.895 0.9478
## block2 - block3 -12.35 47.8 Inf -0.258 0.9998
## block2 - block4 268.10 47.9 Inf 5.602 <.0001
## block2 - block5 270.69 47.7 Inf 5.674 <.0001
## block2 - block6 54.50 48.6 Inf 1.122 0.8724
## block3 - block4 280.45 48.2 Inf 5.820 <.0001
## block3 - block5 283.04 48.0 Inf 5.893 <.0001
## block3 - block6 66.85 48.9 Inf 1.367 0.7467
## block4 - block5 2.59 48.1 Inf 0.054 1.0000
## block4 - block6 -213.59 48.9 Inf -4.365 0.0002
## block5 - block6 -216.19 48.8 Inf -4.430 0.0001
##
## Degrees-of-freedom method: asymptotic
## P value adjustment: tukey method for comparing a family of 6 estimatesae.m.thetaCz<-allEffects(m.thetaCz)
df.ae.m.thetaCz<-as.data.frame(ae.m.thetaCz[[1]])
plot.m.thetaCz <- ggplot(df.ae.m.thetaCz, aes(x=block, y= fit))+
geom_point(aes(color = block)) +
geom_errorbar(aes(ymin= fit-se, ymax= fit+se, fill=block, color=block),
width=.2,position=position_dodge(.9)) +
ylab("ERS/D Theta Power Cz") +
xlab("Block")+
ggtitle("")+
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))## Warning in geom_errorbar(aes(ymin = fit - se, ymax = fit + se, fill = block, :
## Ignoring unknown aesthetics: fillplot(plot.m.thetaCz)
#Theta concat analysis
#factors already made above
m.conthetaCz <- lmer(Cz ~ block * step + (1|subject), data=FINAL_DATA, REML = FALSE)
anova(m.conthetaCz)summary(m.conthetaCz)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: Cz ~ block * step + (1 | subject)
## Data: FINAL_DATA
##
## AIC BIC logLik deviance df.resid
## 302593.9 302887.9 -151258.9 302517.9 16906
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.763 -0.241 -0.075 0.060 36.998
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 449694 670.6
## Residual 3309742 1819.3
## Number of obs: 16944, groups: subject, 12
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 508.237 194.093 12.008 2.619 0.02243 *
## block.L -103.487 35.030 16932.787 -2.954 0.00314 **
## block.Q 56.947 34.879 16932.465 1.633 0.10255
## block.C 241.745 34.145 16932.263 7.080 1.50e-12 ***
## block^4 27.635 33.744 16932.123 0.819 0.41282
## block^5 -88.552 33.784 16932.096 -2.621 0.00877 **
## step.L 29.982 34.287 16932.008 0.874 0.38189
## step.Q -448.154 34.287 16932.008 -13.071 < 2e-16 ***
## step.C -2.902 34.287 16932.008 -0.085 0.93255
## step^4 -14.773 34.287 16932.008 -0.431 0.66657
## step^5 28.892 34.287 16932.008 0.843 0.39943
## block.L:step.L -149.474 85.621 16932.008 -1.746 0.08087 .
## block.Q:step.L 110.439 85.323 16932.008 1.294 0.19556
## block.C:step.L -127.787 83.579 16932.008 -1.529 0.12630
## block^4:step.L 11.048 82.628 16932.008 0.134 0.89363
## block^5:step.L 113.766 82.733 16932.008 1.375 0.16912
## block.L:step.Q 2.753 85.621 16932.008 0.032 0.97435
## block.Q:step.Q -62.292 85.323 16932.008 -0.730 0.46536
## block.C:step.Q -452.047 83.579 16932.008 -5.409 6.44e-08 ***
## block^4:step.Q -41.108 82.628 16932.008 -0.498 0.61884
## block^5:step.Q 114.458 82.733 16932.008 1.383 0.16654
## block.L:step.C 67.710 85.621 16932.008 0.791 0.42906
## block.Q:step.C -71.774 85.323 16932.008 -0.841 0.40025
## block.C:step.C 187.894 83.579 16932.008 2.248 0.02458 *
## block^4:step.C 62.149 82.628 16932.008 0.752 0.45197
## block^5:step.C -79.312 82.733 16932.008 -0.959 0.33775
## block.L:step^4 17.658 85.621 16932.008 0.206 0.83661
## block.Q:step^4 -66.461 85.323 16932.008 -0.779 0.43603
## block.C:step^4 123.848 83.579 16932.008 1.482 0.13841
## block^4:step^4 -24.980 82.628 16932.008 -0.302 0.76241
## block^5:step^4 -62.155 82.733 16932.008 -0.751 0.45250
## block.L:step^5 38.504 85.621 16932.008 0.450 0.65293
## block.Q:step^5 -49.247 85.323 16932.008 -0.577 0.56382
## block.C:step^5 -11.509 83.579 16932.008 -0.138 0.89047
## block^4:step^5 -62.480 82.628 16932.008 -0.756 0.44956
## block^5:step^5 -32.718 82.733 16932.008 -0.395 0.69250
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Correlation matrix not shown by default, as p = 36 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need itemm.conthetaCz <- emmeans(m.conthetaCz, pairwise ~ step|block)## Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
## To enable adjustments, add the argument 'pbkrtest.limit = 16944' (or larger)
## [or, globally, 'set emm_options(pbkrtest.limit = 16944)' 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 = 16944' (or larger)
## [or, globally, 'set emm_options(lmerTest.limit = 16944)' or larger];
## but be warned that this may result in large computation time and memory use.print(emm.conthetaCz)## $emmeans
## block = 1:
## step emmean SE df asymp.LCL asymp.UCL
## 1 241 213 Inf -176.6 658
## 2 461 213 Inf 43.3 878
## 3 575 213 Inf 157.3 992
## 4 675 213 Inf 258.1 1093
## 5 778 213 Inf 361.0 1196
## 6 402 213 Inf -15.7 819
##
## block = 2:
## step emmean SE df asymp.LCL asymp.UCL
## 1 263 210 Inf -149.0 674
## 2 673 210 Inf 261.5 1085
## 3 885 210 Inf 473.6 1297
## 4 932 210 Inf 520.7 1344
## 5 654 210 Inf 242.0 1066
## 6 323 210 Inf -89.1 734
##
## block = 3:
## step emmean SE df asymp.LCL asymp.UCL
## 1 303 211 Inf -109.9 715
## 2 788 211 Inf 375.5 1201
## 3 977 211 Inf 564.6 1390
## 4 912 211 Inf 499.6 1325
## 5 531 211 Inf 118.4 944
## 6 293 211 Inf -119.9 705
##
## block = 4:
## step emmean SE df asymp.LCL asymp.UCL
## 1 185 211 Inf -227.7 598
## 2 370 211 Inf -42.4 783
## 3 396 211 Inf -17.2 808
## 4 459 211 Inf 46.7 872
## 5 465 211 Inf 52.5 878
## 6 246 211 Inf -167.1 659
##
## block = 5:
## step emmean SE df asymp.LCL asymp.UCL
## 1 256 210 Inf -156.2 669
## 2 387 210 Inf -25.4 799
## 3 349 210 Inf -63.8 761
## 4 519 210 Inf 106.4 931
## 5 399 210 Inf -13.7 811
## 6 197 210 Inf -215.8 609
##
## block = 6:
## step emmean SE df asymp.LCL asymp.UCL
## 1 206 212 Inf -208.5 621
## 2 695 212 Inf 279.9 1109
## 3 872 212 Inf 457.4 1287
## 4 827 212 Inf 412.3 1242
## 5 585 212 Inf 170.2 1000
## 6 218 212 Inf -196.3 633
##
## Degrees-of-freedom method: asymptotic
## Confidence level used: 0.95
##
## $contrasts
## block = 1:
## contrast estimate SE df z.ratio p.value
## step1 - step2 -219.91 125 Inf -1.756 0.4946
## step1 - step3 -333.90 125 Inf -2.666 0.0821
## step1 - step4 -434.66 125 Inf -3.471 0.0069
## step1 - step5 -537.56 125 Inf -4.292 0.0003
## step1 - step6 -160.91 125 Inf -1.285 0.7935
## step2 - step3 -113.99 125 Inf -0.910 0.9441
## step2 - step4 -214.75 125 Inf -1.715 0.5219
## step2 - step5 -317.65 125 Inf -2.536 0.1137
## step2 - step6 59.00 125 Inf 0.471 0.9971
## step3 - step4 -100.76 125 Inf -0.804 0.9668
## step3 - step5 -203.65 125 Inf -1.626 0.5812
## step3 - step6 172.99 125 Inf 1.381 0.7384
## step4 - step5 -102.90 125 Inf -0.822 0.9637
## step4 - step6 273.75 125 Inf 2.186 0.2443
## step5 - step6 376.65 125 Inf 3.007 0.0316
##
## block = 2:
## contrast estimate SE df z.ratio p.value
## step1 - step2 -410.55 115 Inf -3.557 0.0050
## step1 - step3 -622.67 115 Inf -5.395 <.0001
## step1 - step4 -669.75 115 Inf -5.803 <.0001
## step1 - step5 -391.04 115 Inf -3.388 0.0092
## step1 - step6 -59.90 115 Inf -0.519 0.9955
## step2 - step3 -212.12 115 Inf -1.838 0.4411
## step2 - step4 -259.20 115 Inf -2.246 0.2168
## step2 - step5 19.51 115 Inf 0.169 1.0000
## step2 - step6 350.65 115 Inf 3.038 0.0287
## step3 - step4 -47.07 115 Inf -0.408 0.9986
## step3 - step5 231.63 115 Inf 2.007 0.3381
## step3 - step6 562.78 115 Inf 4.876 <.0001
## step4 - step5 278.71 115 Inf 2.415 0.1509
## step4 - step6 609.85 115 Inf 5.284 <.0001
## step5 - step6 331.14 115 Inf 2.869 0.0473
##
## block = 3:
## contrast estimate SE df z.ratio p.value
## step1 - step2 -485.40 117 Inf -4.146 0.0005
## step1 - step3 -674.53 117 Inf -5.762 <.0001
## step1 - step4 -609.45 117 Inf -5.206 <.0001
## step1 - step5 -228.33 117 Inf -1.950 0.3713
## step1 - step6 10.01 117 Inf 0.086 1.0000
## step2 - step3 -189.12 117 Inf -1.615 0.5882
## step2 - step4 -124.05 117 Inf -1.060 0.8973
## step2 - step5 257.08 117 Inf 2.196 0.2395
## step2 - step6 495.42 117 Inf 4.232 0.0003
## step3 - step4 65.08 117 Inf 0.556 0.9937
## step3 - step5 446.20 117 Inf 3.811 0.0019
## step3 - step6 684.54 117 Inf 5.847 <.0001
## step4 - step5 381.12 117 Inf 3.256 0.0144
## step4 - step6 619.46 117 Inf 5.291 <.0001
## step5 - step6 238.34 117 Inf 2.036 0.3218
##
## block = 4:
## contrast estimate SE df z.ratio p.value
## step1 - step2 -185.26 117 Inf -1.579 0.6124
## step1 - step3 -210.43 117 Inf -1.794 0.4697
## step1 - step4 -274.34 117 Inf -2.339 0.1786
## step1 - step5 -280.13 117 Inf -2.388 0.1604
## step1 - step6 -60.58 117 Inf -0.516 0.9956
## step2 - step3 -25.17 117 Inf -0.215 0.9999
## step2 - step4 -89.08 117 Inf -0.759 0.9742
## step2 - step5 -94.87 117 Inf -0.809 0.9660
## step2 - step6 124.68 117 Inf 1.063 0.8961
## step3 - step4 -63.91 117 Inf -0.545 0.9943
## step3 - step5 -69.70 117 Inf -0.594 0.9915
## step3 - step6 149.85 117 Inf 1.277 0.7975
## step4 - step5 -5.79 117 Inf -0.049 1.0000
## step4 - step6 213.76 117 Inf 1.822 0.4513
## step5 - step6 219.55 117 Inf 1.872 0.4198
##
## block = 5:
## contrast estimate SE df z.ratio p.value
## step1 - step2 -130.81 117 Inf -1.122 0.8724
## step1 - step3 -92.35 117 Inf -0.792 0.9690
## step1 - step4 -262.60 117 Inf -2.252 0.2140
## step1 - step5 -142.51 117 Inf -1.222 0.8261
## step1 - step6 59.62 117 Inf 0.511 0.9958
## step2 - step3 38.46 117 Inf 0.330 0.9995
## step2 - step4 -131.78 117 Inf -1.130 0.8689
## step2 - step5 -11.70 117 Inf -0.100 1.0000
## step2 - step6 190.43 117 Inf 1.633 0.5762
## step3 - step4 -170.25 117 Inf -1.460 0.6897
## step3 - step5 -50.16 117 Inf -0.430 0.9981
## step3 - step6 151.97 117 Inf 1.304 0.7833
## step4 - step5 120.09 117 Inf 1.030 0.9080
## step4 - step6 322.22 117 Inf 2.764 0.0634
## step5 - step6 202.13 117 Inf 1.734 0.5093
##
## block = 6:
## contrast estimate SE df z.ratio p.value
## step1 - step2 -488.33 121 Inf -4.044 0.0007
## step1 - step3 -665.88 121 Inf -5.515 <.0001
## step1 - step4 -620.76 121 Inf -5.141 <.0001
## step1 - step5 -378.68 121 Inf -3.136 0.0212
## step1 - step6 -12.12 121 Inf -0.100 1.0000
## step2 - step3 -177.54 121 Inf -1.470 0.6833
## step2 - step4 -132.42 121 Inf -1.097 0.8829
## step2 - step5 109.65 121 Inf 0.908 0.9446
## step2 - step6 476.21 121 Inf 3.944 0.0011
## step3 - step4 45.12 121 Inf 0.374 0.9991
## step3 - step5 287.19 121 Inf 2.378 0.1638
## step3 - step6 653.75 121 Inf 5.414 <.0001
## step4 - step5 242.07 121 Inf 2.005 0.3394
## step4 - step6 608.63 121 Inf 5.040 <.0001
## step5 - step6 366.56 121 Inf 3.036 0.0290
##
## Degrees-of-freedom method: asymptotic
## P value adjustment: tukey method for comparing a family of 6 estimatesemmp.m.conthetaCz = emmip(m.conthetaCz, ~ step|block)## Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
## To enable adjustments, add the argument 'pbkrtest.limit = 16944' (or larger)
## [or, globally, 'set emm_options(pbkrtest.limit = 16944)' 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 = 16944' (or larger)
## [or, globally, 'set emm_options(lmerTest.limit = 16944)' or larger];
## but be warned that this may result in large computation time and memory use.print(emmp.m.conthetaCz)
ae.m.conthetaCz<-allEffects(m.conthetaCz)
df.ae.m.conthetaCz<-as.data.frame(ae.m.conthetaCz[[1]])
df.ae.m.conthetaCz$l83 <- df.ae.m.conthetaCz$fit - 1.3722 * df.ae.m.conthetaCz$se
df.ae.m.conthetaCz$u83 <- df.ae.m.conthetaCz$fit + 1.3722 * df.ae.m.conthetaCz$se
plot.m.conthetaCz <- ggplot(df.ae.m.conthetaCz, aes(x=step, y=fit, ymin=l83, ymax=u83))+
geom_point(aes(color = block), size=3) +
geom_path(aes(x=step, y=fit, group=block, colour=block)) +
geom_ribbon(aes(ymin=l83, ymax=u83, group=block, fill=block), alpha=0.2) +
ylab("ERDS (%)")+
xlab("Step Position")+
ggtitle("Concatenation Theta ERDS")+
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))+
facet_wrap(~block)
print(plot.m.conthetaCz)
#Beta frequency #making factors and the model
setwd("C:/Users/Alex/Desktop/eeg data")
eeg_data_beta_RT <- read_excel("C:/Users/Alex/Desktop/eeg data/eeg_data_beta_RT.xlsx")
eeg_data_beta_RT$subject = factor(eeg_data_beta_RT$subject)
eeg_data_beta_RT$block = factor(eeg_data_beta_RT$block, ordered = TRUE, levels=c('1', '2', '3', '4', '5', '6'))
eeg_data_beta_RT$step = factor(eeg_data_beta_RT$step, ordered = TRUE, levels=c('1', '2', '3', '4', '5', '6'))
m.betaCz <- lmer(Cz ~ block + (1|subject), data=eeg_data_beta_RT, REML = FALSE)#AIC BIC
AIC (m.betaCz)## [1] 273579.9BIC (m.betaCz)## [1] 273641.8anova(m.betaCz)summary(m.betaCz)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: Cz ~ block + (1 | subject)
## Data: eeg_data_beta_RT
##
## AIC BIC logLik deviance df.resid
## 273579.9 273641.8 -136781.9 273563.9 16936
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.370 -0.301 -0.138 0.075 44.926
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 55998 236.6
## Residual 599493 774.3
## Number of obs: 16944, groups: subject, 12
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 285.00 68.57 12.01 4.156 0.001329 **
## block.L -38.49 14.98 16933.49 -2.570 0.010185 *
## block.Q -126.02 14.88 16932.64 -8.467 < 2e-16 ***
## block.C -51.07 14.54 16932.44 -3.512 0.000446 ***
## block^4 -73.08 14.34 16932.17 -5.095 3.53e-07 ***
## block^5 -34.66 14.34 16932.07 -2.417 0.015670 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) blck.L blck.Q blck.C blck^4
## block.L -0.002
## block.Q 0.005 -0.048
## block.C -0.002 0.058 -0.026
## block^4 0.002 -0.018 0.022 -0.008
## block^5 0.000 0.013 -0.008 0.000 0.004emm.betaCz <- emmeans(m.betaCz, pairwise ~ block)## Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
## To enable adjustments, add the argument 'pbkrtest.limit = 16944' (or larger)
## [or, globally, 'set emm_options(pbkrtest.limit = 16944)' 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 = 16944' (or larger)
## [or, globally, 'set emm_options(lmerTest.limit = 16944)' or larger];
## but be warned that this may result in large computation time and memory use.print(emm.betaCz)## $emmeans
## block emmean SE df asymp.LCL asymp.UCL
## 1 247 70.1 Inf 109.3 384
## 2 316 69.8 Inf 179.7 453
## 3 324 69.8 Inf 186.8 460
## 4 301 69.8 Inf 164.3 438
## 5 364 69.8 Inf 227.2 501
## 6 158 69.9 Inf 21.3 295
##
## Degrees-of-freedom method: asymptotic
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df z.ratio p.value
## block1 - block2 -69.76 21.2 Inf -3.294 0.0127
## block1 - block3 -76.93 21.3 Inf -3.616 0.0041
## block1 - block4 -54.51 21.4 Inf -2.552 0.1093
## block1 - block5 -117.28 21.3 Inf -5.514 <.0001
## block1 - block6 88.44 21.6 Inf 4.102 0.0006
## block2 - block3 -7.17 20.2 Inf -0.356 0.9993
## block2 - block4 15.25 20.2 Inf 0.754 0.9750
## block2 - block5 -47.52 20.1 Inf -2.360 0.1707
## block2 - block6 158.20 20.4 Inf 7.738 <.0001
## block3 - block4 22.42 20.3 Inf 1.104 0.8800
## block3 - block5 -40.35 20.2 Inf -1.995 0.3449
## block3 - block6 165.37 20.5 Inf 8.057 <.0001
## block4 - block5 -62.77 20.3 Inf -3.093 0.0243
## block4 - block6 142.95 20.6 Inf 6.940 <.0001
## block5 - block6 205.72 20.5 Inf 10.031 <.0001
##
## Degrees-of-freedom method: asymptotic
## P value adjustment: tukey method for comparing a family of 6 estimatesae.m.betaCz<-allEffects(m.betaCz)
df.ae.m.betaCz<-as.data.frame(ae.m.betaCz[[1]])
plot.m.betaCz <- ggplot(df.ae.m.betaCz, aes(x=block, y= fit))+
geom_point(aes(color = block)) +
geom_errorbar(aes(ymin= fit-se, ymax= fit+se, fill=block, color=block),
width=.2,position=position_dodge(.9)) +
ylab("ERS/D Beta Power Cz") +
xlab("Block")+
ggtitle("")+
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))## Warning in geom_errorbar(aes(ymin = fit - se, ymax = fit + se, fill = block, :
## Ignoring unknown aesthetics: fillplot(plot.m.betaCz)
#Beta Concat analysis
#factors already made above
m.conbetaCz <- lmer(Cz ~ block * step + (1|subject), data=eeg_data_beta_RT, REML = FALSE)
anova(m.conbetaCz)summary(m.conbetaCz)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: Cz ~ block * step + (1 | subject)
## Data: eeg_data_beta_RT
##
## AIC BIC logLik deviance df.resid
## 273227.3 273521.4 -136575.7 273151.3 16906
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.606 -0.353 -0.116 0.114 45.300
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 56007 236.7
## Residual 585063 764.9
## Number of obs: 16944, groups: subject, 12
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 285.000 68.573 12.013 4.156 0.001329 **
## block.L -38.490 14.797 16933.460 -2.601 0.009300 **
## block.Q -126.021 14.704 16932.628 -8.571 < 2e-16 ***
## block.C -51.067 14.366 16932.425 -3.555 0.000379 ***
## block^4 -73.077 14.170 16932.165 -5.157 2.53e-07 ***
## block^5 -34.665 14.170 16932.066 -2.446 0.014439 *
## step.L -27.279 14.428 16932.012 -1.891 0.058675 .
## step.Q -238.637 14.428 16932.012 -16.540 < 2e-16 ***
## step.C 74.713 14.428 16932.012 5.178 2.26e-07 ***
## step^4 -59.187 14.428 16932.012 -4.102 4.11e-05 ***
## step^5 29.499 14.428 16932.012 2.045 0.040913 *
## block.L:step.L -16.379 36.143 16932.012 -0.453 0.650424
## block.Q:step.L 106.476 35.973 16932.012 2.960 0.003082 **
## block.C:step.L 62.160 35.160 16932.012 1.768 0.077089 .
## block^4:step.L 66.849 34.697 16932.012 1.927 0.054041 .
## block^5:step.L 65.957 34.705 16932.012 1.901 0.057382 .
## block.L:step.Q 51.632 36.143 16932.012 1.429 0.153151
## block.Q:step.Q 161.067 35.973 16932.012 4.477 7.60e-06 ***
## block.C:step.Q 19.795 35.160 16932.012 0.563 0.573440
## block^4:step.Q 123.009 34.697 16932.012 3.545 0.000393 ***
## block^5:step.Q 25.221 34.705 16932.012 0.727 0.467398
## block.L:step.C -89.988 36.143 16932.012 -2.490 0.012791 *
## block.Q:step.C -74.296 35.973 16932.012 -2.065 0.038908 *
## block.C:step.C 1.647 35.160 16932.012 0.047 0.962639
## block^4:step.C -5.972 34.697 16932.012 -0.172 0.863349
## block^5:step.C -54.748 34.705 16932.012 -1.578 0.114691
## block.L:step^4 44.675 36.143 16932.012 1.236 0.216454
## block.Q:step^4 31.196 35.973 16932.012 0.867 0.385833
## block.C:step^4 -9.537 35.160 16932.012 -0.271 0.786195
## block^4:step^4 -84.060 34.697 16932.012 -2.423 0.015417 *
## block^5:step^4 -6.276 34.705 16932.012 -0.181 0.856490
## block.L:step^5 2.932 36.143 16932.012 0.081 0.935343
## block.Q:step^5 -27.279 35.973 16932.012 -0.758 0.448275
## block.C:step^5 -58.938 35.160 16932.012 -1.676 0.093702 .
## block^4:step^5 -19.924 34.697 16932.012 -0.574 0.565817
## block^5:step^5 -28.214 34.705 16932.012 -0.813 0.416251
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Correlation matrix not shown by default, as p = 36 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need itemm.conbetaCz <- emmeans(m.conbetaCz, pairwise ~ step|block)## Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
## To enable adjustments, add the argument 'pbkrtest.limit = 16944' (or larger)
## [or, globally, 'set emm_options(pbkrtest.limit = 16944)' 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 = 16944' (or larger)
## [or, globally, 'set emm_options(lmerTest.limit = 16944)' or larger];
## but be warned that this may result in large computation time and memory use.print(emm.conbetaCz)## $emmeans
## block = 1:
## step emmean SE df asymp.LCL asymp.UCL
## 1 89.1 78.2 Inf -64.076 242
## 2 358.6 78.2 Inf 205.363 512
## 3 292.1 78.2 Inf 138.900 445
## 4 286.2 78.2 Inf 133.003 439
## 5 262.9 78.2 Inf 109.689 416
## 6 191.1 78.2 Inf 37.863 344
##
## block = 2:
## step emmean SE df asymp.LCL asymp.UCL
## 1 103.7 76.5 Inf -46.189 253
## 2 434.1 76.5 Inf 284.248 584
## 3 475.1 76.5 Inf 325.226 625
## 4 413.8 76.5 Inf 264.002 564
## 5 312.9 76.5 Inf 163.072 463
## 6 159.0 76.5 Inf 9.115 309
##
## block = 3:
## step emmean SE df asymp.LCL asymp.UCL
## 1 134.3 76.6 Inf -15.774 284
## 2 529.8 76.6 Inf 379.696 680
## 3 436.0 76.6 Inf 285.910 586
## 4 372.1 76.6 Inf 221.965 522
## 5 301.5 76.6 Inf 151.414 452
## 6 167.7 76.6 Inf 17.619 318
##
## block = 4:
## step emmean SE df asymp.LCL asymp.UCL
## 1 140.0 76.7 Inf -10.299 290
## 2 436.3 76.7 Inf 285.990 587
## 3 370.2 76.7 Inf 219.859 521
## 4 371.2 76.7 Inf 220.899 522
## 5 333.5 76.7 Inf 183.171 484
## 6 155.6 76.7 Inf 5.295 306
##
## block = 5:
## step emmean SE df asymp.LCL asymp.UCL
## 1 237.4 76.6 Inf 87.385 388
## 2 506.9 76.6 Inf 356.876 657
## 3 493.1 76.6 Inf 343.073 643
## 4 525.1 76.6 Inf 375.068 675
## 5 282.9 76.6 Inf 132.847 433
## 6 138.1 76.6 Inf -11.962 288
##
## block = 6:
## step emmean SE df asymp.LCL asymp.UCL
## 1 76.7 77.0 Inf -74.238 228
## 2 150.1 77.0 Inf -0.844 301
## 3 175.9 77.0 Inf 24.870 327
## 4 190.4 77.0 Inf 39.462 341
## 5 225.5 77.0 Inf 74.484 376
## 6 130.7 77.0 Inf -20.325 282
##
## Degrees-of-freedom method: asymptotic
## Confidence level used: 0.95
##
## $contrasts
## block = 1:
## contrast estimate SE df z.ratio p.value
## step1 - step2 -269.44 53.7 Inf -5.019 <.0001
## step1 - step3 -202.98 53.7 Inf -3.781 0.0022
## step1 - step4 -197.08 53.7 Inf -3.671 0.0033
## step1 - step5 -173.76 53.7 Inf -3.237 0.0153
## step1 - step6 -101.94 53.7 Inf -1.899 0.4028
## step2 - step3 66.46 53.7 Inf 1.238 0.8182
## step2 - step4 72.36 53.7 Inf 1.348 0.7581
## step2 - step5 95.67 53.7 Inf 1.782 0.4773
## step2 - step6 167.50 53.7 Inf 3.120 0.0223
## step3 - step4 5.90 53.7 Inf 0.110 1.0000
## step3 - step5 29.21 53.7 Inf 0.544 0.9943
## step3 - step6 101.04 53.7 Inf 1.882 0.4132
## step4 - step5 23.31 53.7 Inf 0.434 0.9981
## step4 - step6 95.14 53.7 Inf 1.772 0.4839
## step5 - step6 71.83 53.7 Inf 1.338 0.7638
##
## block = 2:
## contrast estimate SE df z.ratio p.value
## step1 - step2 -330.44 48.5 Inf -6.810 <.0001
## step1 - step3 -371.42 48.5 Inf -7.655 <.0001
## step1 - step4 -310.19 48.5 Inf -6.393 <.0001
## step1 - step5 -209.26 48.5 Inf -4.313 0.0002
## step1 - step6 -55.30 48.5 Inf -1.140 0.8648
## step2 - step3 -40.98 48.5 Inf -0.845 0.9591
## step2 - step4 20.25 48.5 Inf 0.417 0.9984
## step2 - step5 121.18 48.5 Inf 2.497 0.1248
## step2 - step6 275.13 48.5 Inf 5.670 <.0001
## step3 - step4 61.22 48.5 Inf 1.262 0.8058
## step3 - step5 162.15 48.5 Inf 3.342 0.0108
## step3 - step6 316.11 48.5 Inf 6.515 <.0001
## step4 - step5 100.93 48.5 Inf 2.080 0.2977
## step4 - step6 254.89 48.5 Inf 5.253 <.0001
## step5 - step6 153.96 48.5 Inf 3.173 0.0189
##
## block = 3:
## contrast estimate SE df z.ratio p.value
## step1 - step2 -395.47 49.0 Inf -8.076 <.0001
## step1 - step3 -301.68 49.0 Inf -6.161 <.0001
## step1 - step4 -237.74 49.0 Inf -4.855 <.0001
## step1 - step5 -167.19 49.0 Inf -3.414 0.0084
## step1 - step6 -33.39 49.0 Inf -0.682 0.9840
## step2 - step3 93.79 49.0 Inf 1.915 0.3926
## step2 - step4 157.73 49.0 Inf 3.221 0.0161
## step2 - step5 228.28 49.0 Inf 4.662 <.0001
## step2 - step6 362.08 49.0 Inf 7.394 <.0001
## step3 - step4 63.95 49.0 Inf 1.306 0.7820
## step3 - step5 134.50 49.0 Inf 2.747 0.0664
## step3 - step6 268.29 49.0 Inf 5.479 <.0001
## step4 - step5 70.55 49.0 Inf 1.441 0.7020
## step4 - step6 204.35 49.0 Inf 4.173 0.0004
## step5 - step6 133.79 49.0 Inf 2.732 0.0690
##
## block = 4:
## contrast estimate SE df z.ratio p.value
## step1 - step2 -296.29 49.3 Inf -6.007 <.0001
## step1 - step3 -230.16 49.3 Inf -4.666 <.0001
## step1 - step4 -231.20 49.3 Inf -4.687 <.0001
## step1 - step5 -193.47 49.3 Inf -3.923 0.0012
## step1 - step6 -15.59 49.3 Inf -0.316 0.9996
## step2 - step3 66.13 49.3 Inf 1.341 0.7622
## step2 - step4 65.09 49.3 Inf 1.320 0.7742
## step2 - step5 102.82 49.3 Inf 2.085 0.2953
## step2 - step6 280.70 49.3 Inf 5.691 <.0001
## step3 - step4 -1.04 49.3 Inf -0.021 1.0000
## step3 - step5 36.69 49.3 Inf 0.744 0.9764
## step3 - step6 214.56 49.3 Inf 4.350 0.0002
## step4 - step5 37.73 49.3 Inf 0.765 0.9733
## step4 - step6 215.60 49.3 Inf 4.371 0.0002
## step5 - step6 177.88 49.3 Inf 3.606 0.0042
##
## block = 5:
## contrast estimate SE df z.ratio p.value
## step1 - step2 -269.49 48.9 Inf -5.515 <.0001
## step1 - step3 -255.69 48.9 Inf -5.232 <.0001
## step1 - step4 -287.68 48.9 Inf -5.887 <.0001
## step1 - step5 -45.46 48.9 Inf -0.930 0.9388
## step1 - step6 99.35 48.9 Inf 2.033 0.3234
## step2 - step3 13.80 48.9 Inf 0.282 0.9998
## step2 - step4 -18.19 48.9 Inf -0.372 0.9991
## step2 - step5 224.03 48.9 Inf 4.584 0.0001
## step2 - step6 368.84 48.9 Inf 7.548 <.0001
## step3 - step4 -32.00 48.9 Inf -0.655 0.9867
## step3 - step5 210.23 48.9 Inf 4.302 0.0002
## step3 - step6 355.03 48.9 Inf 7.265 <.0001
## step4 - step5 242.22 48.9 Inf 4.957 <.0001
## step4 - step6 387.03 48.9 Inf 7.920 <.0001
## step5 - step6 144.81 48.9 Inf 2.963 0.0360
##
## block = 6:
## contrast estimate SE df z.ratio p.value
## step1 - step2 -73.39 50.3 Inf -1.458 0.6909
## step1 - step3 -99.11 50.3 Inf -1.969 0.3601
## step1 - step4 -113.70 50.3 Inf -2.259 0.2110
## step1 - step5 -148.72 50.3 Inf -2.955 0.0369
## step1 - step6 -53.91 50.3 Inf -1.071 0.8929
## step2 - step3 -25.71 50.3 Inf -0.511 0.9958
## step2 - step4 -40.31 50.3 Inf -0.801 0.9674
## step2 - step5 -75.33 50.3 Inf -1.497 0.6664
## step2 - step6 19.48 50.3 Inf 0.387 0.9989
## step3 - step4 -14.59 50.3 Inf -0.290 0.9997
## step3 - step5 -49.61 50.3 Inf -0.986 0.9226
## step3 - step6 45.19 50.3 Inf 0.898 0.9471
## step4 - step5 -35.02 50.3 Inf -0.696 0.9825
## step4 - step6 59.79 50.3 Inf 1.188 0.8428
## step5 - step6 94.81 50.3 Inf 1.884 0.4121
##
## Degrees-of-freedom method: asymptotic
## P value adjustment: tukey method for comparing a family of 6 estimatesae.m.conbetaCz<-allEffects(m.conbetaCz)
df.ae.m.conbetaCz<-as.data.frame(ae.m.conbetaCz[[1]])
df.ae.m.conbetaCz$l83 <- df.ae.m.conbetaCz$fit - 1.3722 * df.ae.m.conbetaCz$se
df.ae.m.conbetaCz$u83 <- df.ae.m.conbetaCz$fit + 1.3722 * df.ae.m.conbetaCz$se
plot.m.conbetaCz <- ggplot(df.ae.m.conbetaCz, aes(x=step, y=fit, ymin=l83, ymax=u83))+
geom_point(aes(color = block), size=3) +
geom_path(aes(x=step, y=fit, group=block, colour=block)) +
geom_ribbon(aes(ymin=l83, ymax=u83, group=block, fill=block), alpha=0.2) +
ylab("ERDS (%)")+
xlab("Step Position")+
ggtitle("Concatenation Beta ERDS")+
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))+
facet_wrap(~block)
print(plot.m.conbetaCz)