Hypotheses
Previous studies showed that the N170 amplitudes evoked by faces preseneted for longer durations were larger:
Duration <- c("17", "200")
Amplitudes <- c(-2.5, -5)
Hypotheses <- rep("means", 2)
df_lit <- tibble(Duration, Amplitudes, Hypotheses)
plot_lit <- {
ggplot(data = df_lit, aes(y = Amplitudes, x = Duration, fill = Duration)) + # set the data, varialbes for x and y axises, and the variable for dividing data
geom_col(color = "black") +
facet_grid(. ~ Hypotheses) +
scale_fill_manual(values = c(gray(.7), gray(.3))) +
scale_y_reverse(limits =c(0, -5.5), breaks = seq(0, -5, -1), expand= c(0, 0)) + # set the limit for y axis
geom_hline(yintercept = -2.5, linetype = "dashed", color = "red", size = 1) +
labs(x = "", y = expression(paste("Amplitudes (", mu, "V)"))) + # set the names for main, x and y axises , title = "Means",
# theme_bw() +
general_theme +
theme( # strip.text = element_text(color = "white"),
legend.position = "none",
axis.line = element_line(size = 3, colour = "black"))
}
plot_lit
There are two possible explainations for the aggregated mean amplitudes:
Duration <- rep(c(rep("17", 5), rep("200", 5) ), 2)
Amplitudes <- c(c(rep(-2.5, 5), rep(-5, 5)), c(rep(-.83, 3), rep(-5, 7)))
TrialNum <- rep(as.character(rep(c(1:5), 2)), 2)
Hypotheses <- c(rep("Possibility 1", 10), rep("Possibility 2", 10))
df_st <- tibble(Duration, Amplitudes, TrialNum, Hypotheses)
plot_st <- {
ggplot(data = df_st, aes(y = Amplitudes, x = Duration, fill = Duration, color = TrialNum)) + # set the data, varialbes for x and y axises, and the variable for dividing data
geom_col(position = "dodge") +
scale_fill_manual(values = c(gray(.7), gray(.3))) +
scale_color_manual(values = rep("black", 5)) +
facet_grid(. ~ Hypotheses) +
geom_hline(yintercept = -2.5, linetype = "dashed", color = "red", size = 1) +
scale_y_reverse(limits =c(0, -5.5), breaks = seq(0, -5, -1), expand= c(0, 0)) + # set the limit for y axis
labs(x = "Duration (ms)") + # set the names for main, x and y axises , title = "At the single-trial level:", , y = expression(paste("Amplitudes (", mu, "V)"))
theme_bw() +
general_theme +
theme(legend.position = "none",
axis.title.y = element_blank(),
axis.line = element_line(size = 3, colour = "black"))
}
plot_st
Brief introduction to the analysis
Equivalence tests
In this study, we do predict that there will be no differences between some conditions. For example, we expect that N170 amplitudes for faces presented for various duration are comparable as long as they are perceived with high subjective confidence. Thus, we employ Equivalence tests to explore if the differences are equivalent to a null region, which we defined as [-0.5, 0.5] in this study. Essentially, the equivalence test performs two one-sided tests to examine if the differences are smaller than the upper boundary of the region and larger than the lower boundary. By performing the conventional null hypothesis significance tests and the equivalence tests, there are (at least) five scenarios for the results:
# equivalence interval
delta_null <- 0.5 # the equivalent interval will be [-delta_null, delta_null]
CI_size = .8
df_equi <- tibble(
Scenarios = paste("Scenario", 1:6),
CI_low = c(delta_null-.15, -delta_null-.5, -delta_null+.15, delta_null-.6, -delta_null-.25, -delta_null-.1),
CI_upp = CI_low + c(rep(CI_size, 5), 1.3)
) %>%
mutate(Scenarios = factor(Scenarios),
Scenarios = factor(Scenarios, levels = rev(levels(Scenarios))))
gg_equi <- ggplot(df_equi, aes(xmin = CI_low, xmax = CI_upp, y = Scenarios)) +
geom_errorbarh(height = .5, color = c("red", "red", "blue", rep("black", 3)), size = 1) +
geom_vline(xintercept = 0, linetype = "longdash") +
geom_vline(xintercept = c(-delta_null, delta_null), linetype = "dashed", color = "gray30") +
scale_x_continuous(limits = c(-1.15, 1.15), breaks = c(-delta_null, 0, delta_null)) +
xlab(expression(paste("Amplitudes (", mu, "V)"))) +
theme_bw() +
theme(
text = element_text(size = 10),
axis.title = element_text(size = 16),
axis.text = element_text(size = 15), # the size of the texts in plot
# axis.text.x = element_text(angle = 45, vjust = 0.5),
legend.title=element_text(size=16),
legend.text=element_text(size=16),
# legend.position = "bottom",
legend.key.width = unit(1.2, "cm"),
plot.title = element_text(lineheight=.8, face="bold", size = 17),
panel.border = element_blank(),
panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
axis.line.x = element_line(colour = 'black', size=0.5, linetype='solid'),
axis.line.y = element_line(colour = 'black', size=0.5, linetype='solid'),
# remove the facet background color
strip.text = element_text(face="bold", size=15, lineheight=5.0),
strip.background = element_rect(fill="white", colour="white", size=1),
strip.placement = "outside"
)
# ggsave("equivalence_test.png", gg_equi, width = 8, height = 4.7)
gg_equi
gg_equi_bw <- ggplot(df_equi, aes(xmin = CI_low, xmax = CI_upp, y = Scenarios)) +
geom_errorbarh(height = .5, size = 1,
color = c("gray50", "gray50", "gray", rep("black", 3))) +
geom_vline(xintercept = 0, linetype = "longdash") +
geom_vline(xintercept = c(-delta_null, delta_null), linetype = "dashed", color = "gray30") +
scale_x_continuous(limits = c(-1.15, 1.15), breaks = c(-delta_null, 0, delta_null)) +
xlab(expression(paste("Amplitudes (", mu, "V)"))) +
theme_bw() +
theme(
text = element_text(size = 10),
axis.title = element_text(size = 16),
axis.text = element_text(size = 15), # the size of the texts in plot
axis.title.y=element_blank(), #remove y axis labels
# axis.text.x = element_text(angle = 45, vjust = 0.5),
legend.title=element_text(size=16),
legend.text=element_text(size=16),
# legend.position = "bottom",
legend.key.width = unit(1.2, "cm"),
plot.title = element_text(lineheight=.8, face="bold", size = 17),
panel.border = element_blank(),
panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
axis.line.x = element_line(colour = 'black', size=0.5, linetype='solid'),
axis.line.y = element_line(colour = 'black', size=0.5, linetype='solid'),
# remove the facet background color
strip.text = element_text(face="bold", size=15, lineheight=5.0),
strip.background = element_rect(fill="white", colour="white", size=1),
strip.placement = "outside"
)
# ggsave("equivalence_test_bw.png", gg_equi_pw, width = 8, height = 4.7)
gg_equi_bw
Behavior results
raw_beha_205 <- {
"E205_35_Behavior.csv" %>%
file.path(folder_data, .) %>%
read_csv()
}
clean_beha_205 <- {
raw_beha_205 %>%
filter(`Procedure[Trial]` == "expProc") %>% # remove rows for practice and countdown
select(ExperimentName, SubjCode = Subject, blockID, Block, Trial, Age, Sex, Handedness, Type = `stimCategory[Trial]`, Category = `FH[Trial]`, Duration = `stimDuration[Trial]`, ACC = `Resp.ACC[Trial]`, Resp.RT = `Resp.RT[Trial]`, Resp = `Resp.RESP[Trial]`, Stimuli = `stimName[Trial]`) %>% # select and rename the columns (remove unuseful columns)
mutate(
SubjCode = factor(paste("P", SubjCode, sep = "")), # save SubjCode as factor
Type = substr(Type, 1, 1),
Type = recode_factor(Type, "N" = "intact", "S" = "scrambled"), # rename the levels of Type
Category = recode(Category, "hosue" = "house", "house" = "house", "face" = "face", .default = "face"), # rename "hosue" as "house"
Duration = factor(Duration), # save Duration as factor
Category = as.factor(Category),
Resp = as.factor(Resp),
blockID = as.factor(blockID)
) %>%
group_by(SubjCode, Type, Category, Duration, ACC) %>% # divide the data based on these conditions
mutate(
RT = Resp.RT + as.numeric(levels(Duration))[Duration],
RT.Z = scale(RT), # calculate the Z value for reaction times within each group
RT.Within3Z = ifelse(RT.Z <= 3 & RT.Z >= -3, 1, ifelse(RT.Z < -3 | RT.Z > 3, 0, NaN)) # if the Z values are between -3 and 3, will be marked as 1
) %>%
ungroup() # ungroup the data
}
# behavior.tidyup
head(clean_beha_205, 10)
# check the number of trials for 205
sum.count.205 <- {
clean_beha_205 %>%
group_by(SubjCode, Type, Category, Duration) %>%
summarize(Count_sum = n())
}
valid.count.205 <- {
clean_beha_205 %>%
filter(RT > RT_min) %>%
group_by(SubjCode, Type, Category, Duration) %>%
summarize(Count = n()) %>%
right_join(sum.count.205, by = c("SubjCode", "Type", "Category", "Duration")) %>%
mutate(ratio = Count / Count_sum) %>%
filter(ratio > ratio_count) %$%
SubjCode %>%
unique()
}
# remove the participants
clean_beha_205 %<>% filter(SubjCode %in% valid.count.205)
# remove participants whose performance in Key 1 was lower than 95%
valid.17 <- {
clean_beha_205 %>%
filter(Duration == 17, Resp == 1) %>% # Type == "intact",
group_by(SubjCode, Duration, Resp) %>% # Type,
summarize(Accuracy = mean(ACC), Count = n()) %>%
filter(Accuracy > 0.95) %$%
SubjCode
}
# clean_beha_205 %>%
# filter(Duration == 17, Resp == 1) %>% # Type == "intact",
# group_by(SubjCode, Duration, Resp) %>% # Type,
# summarize(Accuracy = mean(ACC), Count = n()) %>%
# filter(Accuracy < 0.95) %>%
# ungroup() %>%
# summarize(mean(Accuracy), min(Accuracy), max(Accuracy))
clean_beha_205 %<>% filter(SubjCode %in% valid.17)
Allocation of responses in each condition
tmp.count <- {
clean_beha_205 %>%
group_by(SubjCode, Type, Category, Duration) %>%
summarize(Count_sum = n())
}
avg.dist.long.205 <- {
clean_beha_205 %>%
group_by(SubjCode, Type, Category, Duration, Resp) %>%
summarize(Count = n(), .groups = "drop") %>%
right_join(., tidyr::expand(., SubjCode, Type, Category, Duration, Resp),
by = c("SubjCode", "Type", "Category", "Duration", "Resp")) %>% # make missing rows explicit
mutate(Count = if_else(is.na(Count), 0L, Count)) %>%
right_join(tmp.count, by = c("SubjCode", "Type", "Category", "Duration")) %>%
mutate(RespRate = Count / Count_sum)
}
# descriptive statistics of accuracy for plotting
desc.dist.205 <- {
avg.dist.long.205 %>%
group_by(Type, Category, Duration, Resp) %>%
summarize(N = n()
, Mean = sum(RespRate)/N
, SD = sd(RespRate)
, SE = SD/sqrt(N)
, SE.lo = Mean - SE, SE.hi = Mean + SE
, CI_lo = Mean - SE * 1.96
, CI_hi = Mean + SE * 1.96
, Median = median(RespRate)
, Lower = Mean - SD, Upper = Mean + SD # for rainclound plot
, .groups = "drop"
)
}
avg.dist.wide.205 <- {
avg.dist.long.205 %>%
mutate(Conditions = paste(Type, Category, Duration, Resp, sep = "_")) %>%
dplyr::select(SubjCode, Conditions, RespRate) %>%
spread(Conditions, RespRate)
}
desc.dist.205
dist.RainPlot.205 <- {
ggplot(data = avg.dist.long.205, aes(y = RespRate, x = as.factor(Resp), fill = Category)) +
geom_flat_violin(position = position_nudge(x = .2, y = 0), alpha = .7) +
geom_point(aes(y = RespRate, color = Category), position = position_jitter(width = .15), size = .5, alpha = .8) +
geom_boxplot(aes(y = RespRate), width = .2, outlier.shape = NA, alpha = .7) +
geom_point(data = desc.dist.205, aes(y = Mean, x = Resp, color = Category), position = position_nudge(x = 0.3), size = 2.5) +
geom_errorbar(data = desc.dist.205, aes(y = Mean, ymin = Lower, ymax = Upper), position = position_nudge(x = 0.3), width = 0) +
facet_grid(Type ~ Duration) +
# facet_grid(. ~ Type, switch = "x") + # create two panels to show data
# scale_colour_grey() + # start = .1, end = .6, color for the contour
# scale_fill_grey() + # start = .3, end = .6, color for the fill
# scale_color_brewer(palette = "Set1") + # palette = "RdBu"
# scale_fill_brewer(palette = "Set1") + # palette = "RdBu"
labs(title = "Reaction times for 205", x = "Stimulus Type", y = "Reaction times (ms)", fill = "Duration(ms)", color = "Duration(ms)") + # set the names for main, x and y axises and the legend
# scale_x_discrete(labels=resp_labels) +
# coord_cartesian(ylim = c(0, 1.05)) + # set the limit for y axis
# scale_y_continuous(expand= c(0, 0)) + # remove the space between columns and x axis
theme_bw() + # the backgroud color
plot_theme
}
Â
dura.labels <- c("17ms", "200ms")
names(dura.labels) <- c(17, 200)
dist.ColuPlot.205 <- {
ggplot(data = desc.dist.205, aes(y = Mean, x = Resp, fill = Category)) + # set the data, varialbes for x and y axises, and the variable for dividing data
geom_col(position = "dodge", color = "black", alpha = .7) + # position of columns and countour of columns
facet_grid(Type ~ Duration,
labeller = labeller(Duration = dura.labels)) +
geom_errorbar(mapping = aes(ymin = CI_lo, ymax = CI_hi), linetype = 1, # set the error bar
show.legend = FALSE, width = 0.25, alpha = .5,
position = position_dodge(width=0.9)) +
geom_hline(yintercept = c(1), linetype = 5, alpha = 0.5) + # add the line for 0.5 and 1 (y)
coord_cartesian(ylim = c(0,1.1)) + # set the limit for y axis
scale_y_continuous(expand= c(0, 0)) + # remove the space between columns and x axis
labs(title = "", x = "Responses", y = "Percentage", fill = "Category") + # set the names for main, x and y axises Subjective Responses for E205
scale_fill_grey() + # set the color for columns
# geom_text(label = c("***", "***", "", ""), color = c("red", "red", "red", "red"), size = 6, nudge_y = 0.05, nudge_x = rep(c(-0.225, 0.225), 2)) + # add starts to the significant columns
erp_theme
}
dist.ColuPlot.205
Percentage of face stimuli
df_ratio_205 <- {
clean_beha_205 %>%
mutate(isFace = if_else(Category == "face", 1, 0)) %>%
group_by(SubjCode, Type, Duration, Resp) %>%
summarize(ratio = mean(isFace), Count = n())
}
desc_ratio_205 <- {
df_ratio_205 %>%
group_by(Type, Duration, Resp) %>%
summarize(Mean = mean(ratio)
, N = n()
, SD = sd(ratio)
, SE = SD/sqrt(N)
, SE.lo = Mean - SE, SE.hi = Mean + SE
, CI_lo = Mean - SE * 1.96
, CI_hi = Mean + SE * 1.96
, Median = median(ratio)
, Lower = Mean - SD, Upper = Mean + SD # for rainclound plot
) %>%
ungroup()
}
rate.ColuPlot.205 <- {
ggplot(data = desc_ratio_205, aes(y = Mean, x = Resp)) + # set the data, varialbes for x and y axises, and the variable for dividing data
geom_col(position = "dodge", color = "black", alpha = .7) + # position of columns and countour of columns
facet_grid(Type ~ Duration,
labeller = labeller(Duration = dura.labels)) +
geom_errorbar(mapping = aes(ymin = SE.lo, ymax = SE.hi), linetype = 1, # set the error bar
show.legend = FALSE, width = 0.25, alpha = .5,
position = position_dodge(width=0.9)) +
geom_hline(yintercept = c(0.5, 1), linetype = 5, alpha = 0.5) + # add the line for 0.5 and 1 (y)
# coord_cartesian(ylim = c(0,1.1)) + # set the limit for y axis
scale_y_continuous(expand= c(0, 0), limits = c(-0.1, 1.15), breaks = seq(0, 1, 0.25)) + # remove the space between columns and x axis
labs(x = "Responses", y = "Percentage") + # set the names for main, x and y axises Subjective Responses for E205
scale_fill_grey() + # set the color for columns
# geom_text(label = c("***", "***", "", ""), color = c("red", "red", "red", "red"), size = 6, nudge_y = 0.05, nudge_x = rep(c(-0.225, 0.225), 2)) + # add starts to the significant columns
plot_theme
}
rate.ColuPlot.205
Combine the two plots (for publication):
dura16.labels <- c("33 ms", "216 ms")
names(dura16.labels) <- c(17, 200)
# combine the two behavior response plots
dist_ColuPlot_205 <-
desc.dist.205 %>%
# mutate(Duration = fct_recode(Duration, `33` = "17", `216` = "200")) %>%
ggplot(aes(y = Mean, x = Resp, fill = Category)) + # set the data, varialbes for x and y axises, and the variable for dividing data
geom_col(position = "dodge", color = "black", alpha = .7) + # position of columns and countour of columns
facet_grid(Type ~ Duration,
labeller = labeller(Duration = dura16.labels)) +
geom_errorbar(mapping = aes(ymin = CI_lo, ymax = CI_hi), linetype = 1, # set the error bar
show.legend = FALSE, width = 0.25, alpha = .5,
position = position_dodge(width=0.9)) +
geom_hline(yintercept = c(1), linetype = 5, alpha = 0.5) + # add the line for 0.5 and 1 (y)
# coord_cartesian(ylim = c(0,1.1)) + # set the limit for y axis
scale_y_continuous(expand= c(0, 0), limits = c(-0.1, 1.15), breaks = seq(0, 1, 0.25)) +# remove the space between columns and x axis
labs(x = "Responses", y = "Percentage", fill = "Category") + # set the names for main, x and y axises Subjective Responses for E205
scale_fill_grey() + # set the color for columns
# geom_text(label = c("***", "***", "", ""), color = c("red", "red", "red", "red"), size = 6, nudge_y = 0.05, nudge_x = rep(c(-0.225, 0.225), 2)) + # add starts to the significant columns
erp_theme +
theme(legend.position = c(0.53, 0.3))
resp_ColuPlot_205 <- {
desc_ratio_205 %>%
# mutate(Duration = fct_recode(Duration, "33" = "17", "216" = "200")) %>%
ggplot(aes(y = Mean, x = Resp)) + # set the data, varialbes for x and y axises, and the variable for dividing data
geom_col(position = "dodge", color = "black", alpha = .7) + # position of columns and countour of columns
facet_grid(Type ~ Duration,
labeller = labeller(Duration = dura16.labels)) +
geom_errorbar(mapping = aes(ymin = SE.lo, ymax = SE.hi), linetype = 1, # set the error bar
show.legend = FALSE, width = 0.25, alpha = .5,
position = position_dodge(width=0.9)) +
geom_hline(yintercept = c(0.5, 1), linetype = 5, alpha = 0.5) + # add the line for 0.5 and 1 (y)
# coord_cartesian(ylim = c(0,1.1)) + # set the limit for y axis
scale_y_continuous(expand= c(0, 0), limits = c(-0.1, 1.15), breaks = seq(0, 1, 0.25)) + # remove the space between columns and x axis
labs(x = "Responses", y = "Percentage") + # set the names for main, x and y axises Subjective Responses for E205
scale_fill_grey() + # set the color for columns
# geom_text(label = c("***", "***", "", ""), color = c("red", "red", "red", "red"), size = 6, nudge_y = 0.05, nudge_x = rep(c(-0.225, 0.225), 2)) + # add starts to the significant columns
plot_theme
}
plot_two_response <- ggpubr::ggarrange(dist_ColuPlot_205, resp_ColuPlot_205, nrow = 1,
widths = c(1, .7),
labels = c("A", "B"),
font.label = list(size = 24)) # ratio between the two subplots)
# ggsave('plot_two_response.pdf', plot_two_response, width = 10, height = 5)
ERP amplitude analysis
Load mean amplitude for single trials
The raw mean amplitude of each trial for all participants were loaded.
# read the trial mean amplitude data
df.erp.E2 <- {
"E205_MeanAmp.csv" %>%
file.path(folder_data, .) %>%
read_csv() %>%
substr_Event() %>%
filter(SubjCode %in% valid.17) %>%
mutate(SubjCode = as.factor(SubjCode),
Hemisphere = as.factor(Hemisphere),
urResponse = as.factor(urResponse),
Stimuli = substr(stimName, 2, 6))
}
df.erp.E2 %>%
group_by(SubjCode, Stimuli, Hemisphere, Type, Category, Duration) %>%
summarize(count = n())
Linear mixed model across responses
Set the successive difference coding
# dummy coding
df.erp.E2.all <- sdif_coding_P203_erp(df.erp.E2)
Amplitudes of the P1
# only keep the data for amplitudes of the P1 (correct and incorrect erps)
df.erp.P1.E2 <- {
df.erp.E2.all %>%
filter(Component == "P1")
}
if (saveCSV) {
output.erp.P1.E2 <- file.path(folder_nesi, "E205_erp_P1.RData")
save(df.erp.P1.E2, file = output.erp.P1.E2)
}
# the structure of this dataset
head(df.erp.P1.E2, 10)
The maximal model
file.max.P1.E2 <- file.path(folder_nesi, "E205_P1_lmm_max.RData")
# fit the maximal model
if (!file.exists(file.max.P1.E2)) {
lmm.max.P1.E2 <- lmer(MeanAmp ~ Type * Category * Duration * Hemisphere +
(1 + Type_C + Cate_C + Dura_C + Hemi_C +
Type_Cate + Type_Dura + Cate_Dura + Type_Hemi + Cate_Hemi + Dura_Hemi +
Type_Cate_Dura + Type_Cate_Hemi + Type_Dura_Hemi + Cate_Dura_Hemi +
Type_Cate_Dura_Hemi
| SubjCode) +
(1 + Type_C + Dura_C + Hemi_C +
Type_Dura + Type_Hemi + Dura_Hemi +
Type_Dura_Hemi
| Stimuli),
data = df.erp.P1.E2,
verbose = TRUE,
control = lmerControl(optimizer = "bobyqa", # nloptwrap Nelder_Mead
optCtrl = list(maxfun = 1e7)))
} else {
load(file.max.P1.E2)
}
print(summary(lmm.max.P1.E2), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Type * Category * Duration * Hemisphere + (1 + Type_C + Cate_C + Dura_C + Hemi_C + Type_Cate + Type_Dura + Cate_Dura + Type_Hemi + Cate_Hemi + Dura_Hemi + Type_Cate_Dura + Type_Cate_Hemi + Type_Dura_Hemi + Cate_Dura_Hemi + Type_Cate_Dura_Hemi | SubjCode) + (1 + Type_C + Dura_C + Hemi_C + Type_Dura + Type_Hemi + Dura_Hemi + Type_Dura_Hemi | Stimuli)
## Data: df.erp.P1.E2
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+07))
##
## REML criterion at convergence: 399168.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.4269 -0.5739 -0.0116 0.5710 9.2481
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Stimuli (Intercept) 0.020242 0.14228
## Type_C 0.071284 0.26699 0.05
## Dura_C 0.072281 0.26885 0.40 -0.05
## Hemi_C 0.001384 0.03720 0.70 0.19 0.70
## Type_Dura 0.421342 0.64911 -0.10 0.67 0.12 -0.21
## Type_Hemi 0.011305 0.10632 -0.46 0.71 0.02 0.14 0.32
## Dura_Hemi 0.002092 0.04574 0.49 0.68 0.57 0.49 0.73 0.25
## Type_Dura_Hemi 0.075770 0.27526 -0.34 0.77 0.25 0.03 0.79 0.80 0.63
## SubjCode (Intercept) 1.797320 1.34064
## Type_C 0.113610 0.33706 0.10
## Cate_C 0.050124 0.22388 -0.36 -0.29
## Dura_C 0.321549 0.56705 0.44 -0.09 -0.36
## Hemi_C 0.594791 0.77123 0.31 0.11 0.48 -0.18
## Type_Cate 0.378197 0.61498 0.67 0.25 -0.36 0.38 0.10
## Type_Dura 0.336157 0.57979 0.05 0.22 0.05 0.29 0.02 0.41
## Cate_Dura 0.043316 0.20813 0.54 -0.40 0.09 0.18 0.55 0.41 -0.19
## Type_Hemi 0.034810 0.18657 0.02 0.47 -0.17 -0.02 -0.35 -0.27 -0.25 -0.45
## Cate_Hemi 0.011529 0.10737 0.00 0.55 0.23 -0.23 -0.06 0.03 0.17 -0.61 0.51
## Dura_Hemi 0.239585 0.48947 0.19 -0.08 -0.46 0.46 -0.37 0.23 0.19 0.27 0.22 -0.50
## Type_Cate_Dura 0.382269 0.61828 0.70 0.29 -0.67 0.53 -0.34 0.77 0.23 0.03 0.19 0.22 0.32
## Type_Cate_Hemi 0.127146 0.35658 0.00 -0.38 0.04 0.27 0.07 0.52 0.45 0.46 -0.83 -0.52 0.20 0.06
## Type_Dura_Hemi 0.133551 0.36545 -0.48 0.43 0.55 -0.36 0.34 -0.39 0.23 -0.26 0.30 0.29 -0.07 -0.59 -0.26
## Cate_Dura_Hemi 0.074409 0.27278 0.06 -0.44 0.18 0.21 -0.41 -0.16 -0.04 0.05 0.46 0.04 0.50 0.10 -0.18 -0.03
## Type_Cate_Dura_Hemi 0.470022 0.68558 0.17 -0.64 0.44 0.34 -0.05 0.17 0.28 0.38 -0.16 -0.10 0.30 0.07 0.38 -0.12 0.75
## Residual 15.678168 3.95957
## Number of obs: 71282, groups: Stimuli, 80; SubjCode, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.262713 0.245906 29.244440 5.135 1.7e-05 ***
## Type2-1 -0.106661 0.076825 37.967608 -1.388 0.17312
## Category2-1 0.085615 0.062504 46.064875 1.370 0.17741
## Duration2-1 0.353467 0.113340 33.150333 3.119 0.00374 **
## Hemisphere2-1 0.117784 0.145144 29.038601 0.811 0.42368
## Type2-1:Category2-1 0.122340 0.145142 38.767656 0.843 0.40445
## Type2-1:Duration2-1 0.313789 0.146180 44.709815 2.147 0.03729 *
## Category2-1:Duration2-1 0.122743 0.099766 60.821395 1.230 0.22332
## Type2-1:Hemisphere2-1 0.011121 0.078696 72.888288 0.141 0.88801
## Category2-1:Hemisphere2-1 -0.041851 0.073108 178.077600 -0.572 0.56774
## Duration2-1:Hemisphere2-1 0.005918 0.113598 30.243456 0.052 0.95879
## Type2-1:Category2-1:Duration2-1 0.008769 0.231065 63.040636 0.038 0.96985
## Type2-1:Category2-1:Hemisphere2-1 0.190520 0.156105 79.923915 1.220 0.22588
## Type2-1:Duration2-1:Hemisphere2-1 -0.044562 0.158002 76.945281 -0.282 0.77867
## Category2-1:Duration2-1:Hemisphere2-1 -0.151587 0.148829 122.513631 -1.019 0.31043
## Type2-1:Category2-1:Duration2-1:Hemisphere2-1 0.123037 0.312596 84.251373 0.394 0.69487
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## optimizer (bobyqa) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
The zero-correlation-parameter model
file.zcp.P1.E2 <- file.path(folder_nesi, "E205_P1_lmm_zcp.RData")
# fit the zcp model
if (!file.exists(file.zcp.P1.E2)) {
lmm.zcp.P1.E2 <- lmer(MeanAmp ~ Type * Category * Duration * Hemisphere +
(1 + Type_C + Cate_C + Dura_C + Hemi_C +
Type_Cate + Type_Dura + Cate_Dura + Type_Hemi + Cate_Hemi + Dura_Hemi +
Type_Cate_Dura + Type_Cate_Hemi + Type_Dura_Hemi + Cate_Dura_Hemi +
Type_Cate_Dura_Hemi
|| SubjCode) +
(1 + Type_C + Dura_C + Hemi_C +
Type_Dura + Type_Hemi + Dura_Hemi +
Type_Dura_Hemi
|| Stimuli),
data = df.erp.P1.E2,
verbose = TRUE,
control = lmerControl(optimizer = "bobyqa", # nloptwrap Nelder_Mead
optCtrl = list(maxfun = 1e7)))
} else {
load(file.zcp.P1.E2)
}
print(summary(lmm.zcp.P1.E2), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Type * Category * Duration * Hemisphere + (1 + Type_C + Cate_C + Dura_C + Hemi_C + Type_Cate + Type_Dura + Cate_Dura + Type_Hemi + Cate_Hemi + Dura_Hemi + Type_Cate_Dura + Type_Cate_Hemi + Type_Dura_Hemi + Cate_Dura_Hemi + Type_Cate_Dura_Hemi || SubjCode) + (1 + Type_C + Dura_C + Hemi_C + Type_Dura + Type_Hemi + Dura_Hemi + Type_Dura_Hemi || Stimuli)
## Data: df.erp.P1.E2
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+07))
##
## REML criterion at convergence: 399290.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.5143 -0.5726 -0.0123 0.5716 9.3014
##
## Random effects:
## Groups Name Variance Std.Dev.
## Stimuli Type_Dura_Hemi 1.932e-09 4.395e-05
## Stimuli.1 Dura_Hemi 2.685e-13 5.181e-07
## Stimuli.2 Type_Hemi 4.559e-13 6.752e-07
## Stimuli.3 Type_Dura 2.496e-01 4.996e-01
## Stimuli.4 Hemi_C 0.000e+00 0.000e+00
## Stimuli.5 Dura_C 5.021e-02 2.241e-01
## Stimuli.6 Type_C 3.631e-02 1.906e-01
## Stimuli.7 (Intercept) 1.472e-02 1.213e-01
## SubjCode Type_Cate_Dura_Hemi 3.519e-12 1.876e-06
## SubjCode.1 Cate_Dura_Hemi 1.916e-13 4.377e-07
## SubjCode.2 Type_Dura_Hemi 0.000e+00 0.000e+00
## SubjCode.3 Type_Cate_Hemi 0.000e+00 0.000e+00
## SubjCode.4 Type_Cate_Dura 0.000e+00 0.000e+00
## SubjCode.5 Dura_Hemi 2.525e-01 5.025e-01
## SubjCode.6 Cate_Hemi 0.000e+00 0.000e+00
## SubjCode.7 Type_Hemi 0.000e+00 0.000e+00
## SubjCode.8 Cate_Dura 3.003e-02 1.733e-01
## SubjCode.9 Type_Dura 3.098e-01 5.566e-01
## SubjCode.10 Type_Cate 2.762e-01 5.255e-01
## SubjCode.11 Hemi_C 6.411e-01 8.007e-01
## SubjCode.12 Dura_C 3.127e-01 5.592e-01
## SubjCode.13 Cate_C 4.097e-02 2.024e-01
## SubjCode.14 Type_C 1.055e-01 3.248e-01
## SubjCode.15 (Intercept) 1.781e+00 1.334e+00
## Residual 1.570e+01 3.963e+00
## Number of obs: 71282, groups: Stimuli, 80; SubjCode, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.263e+00 2.446e-01 2.918e+01 5.162 1.59e-05 ***
## Type2-1 -1.068e-01 7.208e-02 3.472e+01 -1.482 0.1474
## Category2-1 8.635e-02 5.768e-02 4.827e+01 1.497 0.1409
## Duration2-1 3.532e-01 1.108e-01 3.198e+01 3.187 0.0032 **
## Hemisphere2-1 1.172e-01 1.503e-01 2.930e+01 0.780 0.4419
## Type2-1:Category2-1 1.227e-01 1.262e-01 4.055e+01 0.972 0.3366
## Type2-1:Duration2-1 3.127e-01 1.355e-01 4.016e+01 2.309 0.0262 *
## Category2-1:Duration2-1 1.223e-01 9.173e-02 4.585e+01 1.333 0.1891
## Type2-1:Hemisphere2-1 1.154e-02 6.999e-02 7.068e+04 0.165 0.8691
## Category2-1:Hemisphere2-1 -4.049e-02 6.999e-02 7.068e+04 -0.579 0.5629
## Duration2-1:Hemisphere2-1 5.761e-03 1.154e-01 3.009e+01 0.050 0.9605
## Type2-1:Category2-1:Duration2-1 7.014e-03 1.791e-01 9.725e+01 0.039 0.9688
## Type2-1:Category2-1:Hemisphere2-1 1.898e-01 1.400e-01 7.068e+04 1.356 0.1752
## Type2-1:Duration2-1:Hemisphere2-1 -4.251e-02 1.400e-01 7.068e+04 -0.304 0.7613
## Category2-1:Duration2-1:Hemisphere2-1 -1.505e-01 1.400e-01 7.068e+04 -1.075 0.2822
## Type2-1:Category2-1:Duration2-1:Hemisphere2-1 1.214e-01 2.800e-01 7.068e+04 0.434 0.6646
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## optimizer (bobyqa) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
The reduced model
summary(rePCA(lmm.zcp.P1.E2))
## $Stimuli
## Importance of components:
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## Standard deviation 0.1261 0.05655 0.04809 0.03062 1.109e-05 1.704e-07 1.308e-07 0
## Proportion of Variance 0.7114 0.14312 0.10350 0.04196 0.000e+00 0.000e+00 0.000e+00 0
## Cumulative Proportion 0.7114 0.85454 0.95804 1.00000 1.000e+00 1.000e+00 1.000e+00 1
##
## $SubjCode
## Importance of components:
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16]
## Standard deviation 0.3368 0.2021 0.14111 0.14047 0.13262 0.12681 0.08195 0.05108 0.04373 4.734e-07 1.105e-07 0 0 0 0 0
## Proportion of Variance 0.4749 0.1710 0.08339 0.08264 0.07366 0.06734 0.02813 0.01093 0.00801 0.000e+00 0.000e+00 0 0 0 0 0
## Cumulative Proportion 0.4749 0.6459 0.72930 0.81194 0.88559 0.95294 0.98106 0.99199 1.00000 1.000e+00 1.000e+00 1 1 1 1 1
file.rdc.P1.E2 <- file.path(folder_nesi, "E205_P1_lmm_rdc.RData")
# fit the rdc model
if (!file.exists(file.rdc.P1.E2)) {
# lmm.rdc.P1.E2.step <- step(lmm.rdc.P1.E2, reduce.fixed = FALSE)
lmm.rdc.P1.E2 <- update(lmm.zcp.P1.E2,
formula = MeanAmp ~ Type * Category * Duration * Hemisphere +
(1 + Type_C + Cate_C + Dura_C + Hemi_C +
Type_Cate + Type_Dura + Cate_Dura + Dura_Hemi
|| SubjCode) +
(1 + Type_C + Dura_C +
Type_Dura
|| Stimuli))
} else {
load(file.rdc.P1.E2)
}
print(summary(lmm.rdc.P1.E2), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Type + Category + Duration + Hemisphere + (1 + Type_C + Cate_C + Dura_C + Hemi_C + Type_Cate + Type_Dura + Cate_Dura + Dura_Hemi || SubjCode) + (1 + Type_C + Dura_C + Type_Dura || Stimuli) + Type:Category + Type:Duration + Category:Duration + Type:Hemisphere + Category:Hemisphere + Duration:Hemisphere + Type:Category:Duration + Type:Category:Hemisphere + Type:Duration:Hemisphere + Category:Duration:Hemisphere + Type:Category:Duration:Hemisphere
## Data: df.erp.P1.E2
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+07))
##
## REML criterion at convergence: 399290.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.5143 -0.5726 -0.0123 0.5716 9.3014
##
## Random effects:
## Groups Name Variance Std.Dev.
## Stimuli Type_Dura 0.24958 0.4996
## Stimuli.1 Dura_C 0.05021 0.2241
## Stimuli.2 Type_C 0.03631 0.1906
## Stimuli.3 (Intercept) 0.01472 0.1213
## SubjCode Dura_Hemi 0.25250 0.5025
## SubjCode.1 Cate_Dura 0.03003 0.1733
## SubjCode.2 Type_Dura 0.30984 0.5566
## SubjCode.3 Type_Cate 0.27617 0.5255
## SubjCode.4 Hemi_C 0.64112 0.8007
## SubjCode.5 Dura_C 0.31265 0.5592
## SubjCode.6 Cate_C 0.04097 0.2024
## SubjCode.7 Type_C 0.10546 0.3248
## SubjCode.8 (Intercept) 1.78060 1.3344
## Residual 15.70241 3.9626
## Number of obs: 71282, groups: Stimuli, 80; SubjCode, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.263e+00 2.446e-01 2.918e+01 5.162 1.59e-05 ***
## Type2-1 -1.068e-01 7.208e-02 3.472e+01 -1.482 0.1474
## Category2-1 8.635e-02 5.768e-02 4.827e+01 1.497 0.1409
## Duration2-1 3.532e-01 1.108e-01 3.198e+01 3.187 0.0032 **
## Hemisphere2-1 1.172e-01 1.503e-01 2.930e+01 0.780 0.4419
## Type2-1:Category2-1 1.227e-01 1.262e-01 4.055e+01 0.972 0.3366
## Type2-1:Duration2-1 3.127e-01 1.355e-01 4.016e+01 2.309 0.0262 *
## Category2-1:Duration2-1 1.223e-01 9.173e-02 4.585e+01 1.333 0.1891
## Type2-1:Hemisphere2-1 1.154e-02 6.999e-02 7.068e+04 0.165 0.8691
## Category2-1:Hemisphere2-1 -4.049e-02 6.999e-02 7.068e+04 -0.579 0.5629
## Duration2-1:Hemisphere2-1 5.761e-03 1.154e-01 3.009e+01 0.050 0.9605
## Type2-1:Category2-1:Duration2-1 7.014e-03 1.791e-01 9.725e+01 0.039 0.9688
## Type2-1:Category2-1:Hemisphere2-1 1.898e-01 1.400e-01 7.068e+04 1.356 0.1752
## Type2-1:Duration2-1:Hemisphere2-1 -4.251e-02 1.400e-01 7.068e+04 -0.304 0.7613
## Category2-1:Duration2-1:Hemisphere2-1 -1.505e-01 1.400e-01 7.068e+04 -1.075 0.2822
## Type2-1:Category2-1:Duration2-1:Hemisphere2-1 1.214e-01 2.800e-01 7.068e+04 0.434 0.6646
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lmm.rdc.P1.E2, lmm.zcp.P1.E2, refit = FALSE)
The extended model
file.etd.P1.E2 <- file.path(folder_nesi, "E205_P1_lmm_etd.RData")
# fit the etd model
if (!file.exists(file.etd.P1.E2)) {
lmm.etd.P1.E2 <- update(lmm.rdc.P1.E2,
formula = MeanAmp ~ Type * Category * Duration * Hemisphere +
(1 + Type_C + Cate_C + Dura_C + Hemi_C +
Type_Cate + Type_Dura + Cate_Dura + Dura_Hemi
| SubjCode) +
(1 + Type_C + Dura_C +
Type_Dura
| Stimuli))
lmm.etd.P1.E2.2 <- re_fit(lmm.etd.P1.E2)
} else {
load(file.etd.P1.E2)
}
print(summary(lmm.etd.P1.E2.2), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Type + Category + Duration + Hemisphere + (1 + Type_C + Cate_C + Dura_C + Hemi_C + Type_Cate + Type_Dura + Cate_Dura + Dura_Hemi | SubjCode) + (1 + Type_C + Dura_C + Type_Dura | Stimuli) + Type:Category + Type:Duration + Category:Duration + Type:Hemisphere + Category:Hemisphere + Duration:Hemisphere + Type:Category:Duration + Type:Category:Hemisphere + Type:Duration:Hemisphere + Category:Duration:Hemisphere + Type:Category:Duration:Hemisphere
## Data: df.erp.P1.E2
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+07))
##
## REML criterion at convergence: 399219.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.4381 -0.5727 -0.0121 0.5713 9.2907
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Stimuli (Intercept) 0.02013 0.1419
## Type_C 0.06908 0.2628 0.06
## Dura_C 0.07107 0.2666 0.40 -0.06
## Type_Dura 0.41428 0.6436 -0.10 0.67 0.12
## SubjCode (Intercept) 1.79762 1.3408
## Type_C 0.11245 0.3353 0.11
## Cate_C 0.04079 0.2020 -0.35 -0.31
## Dura_C 0.32097 0.5665 0.44 -0.09 -0.37
## Hemi_C 0.61141 0.7819 0.30 0.11 0.50 -0.19
## Type_Cate 0.25555 0.5055 0.62 0.23 -0.23 0.31 0.22
## Type_Dura 0.33290 0.5770 0.05 0.22 0.07 0.29 0.03 0.43
## Cate_Dura 0.04222 0.2055 0.45 -0.45 0.13 0.10 0.58 0.41 -0.23
## Dura_Hemi 0.23709 0.4869 0.22 -0.08 -0.60 0.48 -0.43 0.20 0.16 0.14
## Residual 15.69460 3.9616
## Number of obs: 71282, groups: Stimuli, 80; SubjCode, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.263e+00 2.459e-01 2.924e+01 5.135 1.7e-05 ***
## Type2-1 -1.065e-01 7.640e-02 3.752e+01 -1.394 0.17163
## Category2-1 8.568e-02 5.993e-02 4.572e+01 1.430 0.15958
## Duration2-1 3.540e-01 1.132e-01 3.306e+01 3.127 0.00367 **
## Hemisphere2-1 1.175e-01 1.470e-01 2.908e+01 0.799 0.43066
## Type2-1:Category2-1 1.216e-01 1.299e-01 4.603e+01 0.936 0.35391
## Type2-1:Duration2-1 3.135e-01 1.455e-01 4.433e+01 2.154 0.03669 *
## Category2-1:Duration2-1 1.222e-01 9.930e-02 6.054e+01 1.231 0.22312
## Type2-1:Hemisphere2-1 1.162e-02 6.997e-02 7.073e+04 0.166 0.86805
## Category2-1:Hemisphere2-1 -4.073e-02 6.997e-02 7.073e+04 -0.582 0.56048
## Duration2-1:Hemisphere2-1 6.368e-03 1.131e-01 2.998e+01 0.056 0.95549
## Type2-1:Category2-1:Duration2-1 6.142e-03 2.008e-01 7.769e+01 0.031 0.97567
## Type2-1:Category2-1:Hemisphere2-1 1.886e-01 1.399e-01 7.073e+04 1.348 0.17771
## Type2-1:Duration2-1:Hemisphere2-1 -4.223e-02 1.399e-01 7.073e+04 -0.302 0.76282
## Category2-1:Duration2-1:Hemisphere2-1 -1.508e-01 1.399e-01 7.073e+04 -1.078 0.28117
## Type2-1:Category2-1:Duration2-1:Hemisphere2-1 1.195e-01 2.799e-01 7.073e+04 0.427 0.66937
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lmm.etd.P1.E2.2, lmm.rdc.P1.E2, refit = FALSE)
The reduced model (lmm.rdc.P1.E2) explained the data better than the extended model (lmm.rdc.P1.E2.2).
The optimal model
lmm.opt.P1.E205 <- lmm.rdc.P1.E2
# print(summary(lmm.opt.P1.E2), corr = FALSE)
summary(lmm.opt.P1.E205)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Type + Category + Duration + Hemisphere + (1 + Type_C + Cate_C + Dura_C + Hemi_C + Type_Cate + Type_Dura + Cate_Dura + Dura_Hemi || SubjCode) + (1 + Type_C + Dura_C + Type_Dura || Stimuli) + Type:Category + Type:Duration + Category:Duration + Type:Hemisphere + Category:Hemisphere + Duration:Hemisphere + Type:Category:Duration + Type:Category:Hemisphere + Type:Duration:Hemisphere + Category:Duration:Hemisphere + Type:Category:Duration:Hemisphere
## Data: df.erp.P1.E2
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+07))
##
## REML criterion at convergence: 399290.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.5143 -0.5726 -0.0123 0.5716 9.3014
##
## Random effects:
## Groups Name Variance Std.Dev.
## Stimuli Type_Dura 0.24958 0.4996
## Stimuli.1 Dura_C 0.05021 0.2241
## Stimuli.2 Type_C 0.03631 0.1906
## Stimuli.3 (Intercept) 0.01472 0.1213
## SubjCode Dura_Hemi 0.25250 0.5025
## SubjCode.1 Cate_Dura 0.03003 0.1733
## SubjCode.2 Type_Dura 0.30984 0.5566
## SubjCode.3 Type_Cate 0.27617 0.5255
## SubjCode.4 Hemi_C 0.64112 0.8007
## SubjCode.5 Dura_C 0.31265 0.5592
## SubjCode.6 Cate_C 0.04097 0.2024
## SubjCode.7 Type_C 0.10546 0.3248
## SubjCode.8 (Intercept) 1.78060 1.3344
## Residual 15.70241 3.9626
## Number of obs: 71282, groups: Stimuli, 80; SubjCode, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.263e+00 2.446e-01 2.918e+01 5.162 1.59e-05 ***
## Type2-1 -1.068e-01 7.208e-02 3.472e+01 -1.482 0.1474
## Category2-1 8.635e-02 5.768e-02 4.827e+01 1.497 0.1409
## Duration2-1 3.532e-01 1.108e-01 3.198e+01 3.187 0.0032 **
## Hemisphere2-1 1.172e-01 1.503e-01 2.930e+01 0.780 0.4419
## Type2-1:Category2-1 1.227e-01 1.262e-01 4.055e+01 0.972 0.3366
## Type2-1:Duration2-1 3.127e-01 1.355e-01 4.016e+01 2.309 0.0262 *
## Category2-1:Duration2-1 1.223e-01 9.173e-02 4.585e+01 1.333 0.1891
## Type2-1:Hemisphere2-1 1.154e-02 6.999e-02 7.068e+04 0.165 0.8691
## Category2-1:Hemisphere2-1 -4.049e-02 6.999e-02 7.068e+04 -0.579 0.5629
## Duration2-1:Hemisphere2-1 5.761e-03 1.154e-01 3.009e+01 0.050 0.9605
## Type2-1:Category2-1:Duration2-1 7.014e-03 1.791e-01 9.725e+01 0.039 0.9688
## Type2-1:Category2-1:Hemisphere2-1 1.898e-01 1.400e-01 7.068e+04 1.356 0.1752
## Type2-1:Duration2-1:Hemisphere2-1 -4.251e-02 1.400e-01 7.068e+04 -0.304 0.7613
## Category2-1:Duration2-1:Hemisphere2-1 -1.505e-01 1.400e-01 7.068e+04 -1.075 0.2822
## Type2-1:Category2-1:Duration2-1:Hemisphere2-1 1.214e-01 2.800e-01 7.068e+04 0.434 0.6646
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# confint(lmm.opt.P1.E205, method="Wald")
Diagnostic plots
# qqplot
qqplot_lmer(lmm.opt.P1.E205)
Estimated marginal means
emm.P1.E205 <- emmeans(lmm.opt.P1.E205, ~ Type + Hemisphere + Category + Duration)
emm.P1.E205 %>%
as.data.frame() %>%
arrange(Type, Hemisphere, Category, Duration)
Plots
# emmip(emm.P1.E205, Category ~ Duration | Type + Hemisphere, CIs = TRUE)
hemiLabel <- c("left", "right")
names(hemiLabel) <- c("Left", "Right")
p1.all.LinePlot = {
ggplot(data = data.frame(emm.P1.E205), aes(y = emmean, x = Duration, color = Category, group = Category)) + # set the data, varialbes for x and y axises, and the variable for dividing data
geom_point(size = 2) +
geom_line(aes(linetype = Category), size = .85) + # position = "dodge", alpha = .7
scale_linetype_manual(values=c("solid", "dashed")) +
geom_errorbar(mapping = aes(ymin = lower.CL, ymax = upper.CL), linetype = 1, # set the error bar
show.legend = FALSE, width = 0.25, alpha = .5) + # , position = position_dodge(width=0.9)
facet_grid(Type ~ Hemisphere, # switch = "x",
labeller = labeller(Hemisphere = hemiLabel)) +
scale_y_reverse(limits = ylimit_p1) + # set the limit for y axis
# scale_y_continuous(expand= c(0, 0)) + # remove the space between columns and x axis
labs(x = "Duration (ms)", y = expression(paste("P1 Amplitudes (", mu, "V)")), fill = "Category") + # set the names for main, x and y axises title = "",
# scale_x_discrete(labels = c("left", "right")) +
# geom_text(label = c("", "", "", "", "***", "", "", ""), color = "red", size = 6, nudge_y = 0.5, nudge_x = 0.5) + # add starts to the significant columns
theme_bw() +
erp_theme
}
# p1.all.LinePlot.slide <- {
# p1.all.LinePlot +
# theme(legend.position="right",
# axis.line = element_line(color="black", size = 3))
# }
# ggsave('P1_all.png', p1.all.LinePlot.slide, width = 7, height = 4)
p1.all.LinePlot
library(RColorBrewer)
duration_color <- c(brewer.pal(9, 'Blues')[c(6:8)], brewer.pal(12, "Paired")[6])
# category_color <- brewer.pal(8, 'Dark2')[c(6,1)]
p1_all_LinePlot <- {
data.frame(emm.P1.E205) %>%
mutate(Duration = fct_recode(Duration, `33`="17", `216`="200")) %>%
ggplot(aes(y = emmean, x = Duration, group = Category)) + # set the data, varialbes for x and y axises, and the variable for dividing data
geom_point(aes(fill = Duration), size=4, shape=21, stroke=0) +
scale_fill_manual(values = duration_color[c(3, 4)], guide = "none") +
geom_line(aes(linetype = Category)) + # position = "dodge", alpha = .7
# scale_color_manual(values=category_color, limits = c("face", "house")) +
scale_linetype_manual(values=c("solid", "dashed")) +
geom_errorbar(mapping = aes(ymin = lower.CL, ymax = upper.CL), linetype = 1, # set the error bar
show.legend = FALSE, width = 0.25, alpha = .5) + # , position = position_dodge(width=0.9)
facet_grid(Type ~ Hemisphere, # switch = "x",
labeller = labeller(Hemisphere = hemiLabel)) +
scale_y_reverse(limits = ylimit_p1) + # set the limit for y axis
# scale_y_continuous(expand= c(0, 0)) + # remove the space between columns and x axis
labs(x = "Duration (ms)", y = expression(paste("P1 Amplitudes (", mu, "V)")), color = "Category", linetype = "Category") + # set the names for main, x and y axises title = "",
# scale_x_discrete(labels = c("left", "right")) +
# geom_text(label = c("", "", "", "", "***", "", "", ""), color = "red", size = 6, nudge_y = 0.5, nudge_x = 0.5) + # add starts to the significant columns
theme_bw() +
erp_theme +
theme(
legend.position = c(.5, .5),
legend.title=element_text(size=10),
legend.text=element_text(size=10),
) +
guides(linetype = guide_legend(title.position = "left", nrow = 1),
color = guide_legend(title.position = "left", nrow = 1))
}
ggsave('P1_all_pub.pdf', p1_all_LinePlot, width = 6, height = 5)
Contrasts
Main results for intact stimuli
Comparisons of P1 amplitudes for intact faces with different durations:
contra_all_P1 <- contrast(emmeans(lmm.opt.P1.E205, ~ Duration | Hemisphere + Category + Type),
"pairwise", combine = TRUE, adjust = "Bonferroni", infer = TRUE)
contra_all_P1[1:2]
## contrast Hemisphere Category Type estimate SE df lower.CL upper.CL t.ratio p.value
## 17 - 200 Left face intact -0.07107583 0.1643558 122.97 -0.4440288 0.3018771 -0.432 1.0000
## 17 - 200 Right face intact -0.20370672 0.1643558 122.97 -0.5766597 0.1692462 -1.239 0.4351
##
## Degrees-of-freedom method: satterthwaite
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 2 estimates
## P value adjustment: bonferroni method for 2 tests
summary(contra_all_P1[1:2], delta = delta_null, level = .9)
summary(contra_all_P1[1:2], delta = delta_null, side = ">")
summary(contra_all_P1[1:2], delta = delta_null, side = "<")
Comparisons of face-specific P1 amplitude components for intact stimuli with different durations:
contra_all_P1_spec <- contrast(emmeans(lmm.opt.P1.E205, ~ Category + Duration | Hemisphere + Type),
interaction = "pairwise", adjust = "Bonferroni", infer = TRUE)
contra_all_P1_spec[1:2]
## Category_pairwise Duration_pairwise Hemisphere Type estimate SE df lower.CL upper.CL t.ratio p.value
## face - house 17 - 200 Left intact 0.22438787 0.1561692 266.23 -0.1276411 0.5764168 1.437 0.3039
## face - house 17 - 200 Right intact 0.01316323 0.1561692 266.23 -0.3388657 0.3651922 0.084 1.0000
##
## Degrees-of-freedom method: satterthwaite
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 2 estimates
## P value adjustment: bonferroni method for 2 tests
summary(contra_all_P1_spec[1:2], delta = delta_null, level = .9)
summary(contra_all_P1_spec[1:2], delta = delta_null, side = ">")
summary(contra_all_P1_spec[1:2], delta = delta_null, side = "<")
Main results for scrambled stimuli
Comparisons of P1 amplitudes for scrambled faces with different durations:
contra_all_P1[5:6]
## contrast Hemisphere Category Type estimate SE df lower.CL upper.CL t.ratio p.value
## 17 - 200 Left face scrambled -0.4319184 0.1700147 140.79 -0.8171081 -0.04672872 -2.540 0.0243
## 17 - 200 Right face scrambled -0.4613435 0.1700147 140.79 -0.8465331 -0.07615380 -2.714 0.0150
##
## Degrees-of-freedom method: satterthwaite
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 2 estimates
## P value adjustment: bonferroni method for 2 tests
summary(contra_all_P1[5:6], delta = delta_null, level = .9)
summary(contra_all_P1[5:6], delta = delta_null, side = ">")
summary(contra_all_P1[5:6], delta = delta_null, side = "<")
Comparisons of face-specific P1 amplitude components for scrambled stimuli with different durations:
contra_all_P1_spec[3:4]
## Category_pairwise Duration_pairwise Hemisphere Type estimate SE df lower.CL upper.CL t.ratio p.value
## face - house 17 - 200 Left scrambled 0.17071189 0.1675616 352.74 -0.2064712 0.5478950 1.019 0.6180
## face - house 17 - 200 Right scrambled 0.08086813 0.1675616 352.74 -0.2963150 0.4580513 0.483 1.0000
##
## Degrees-of-freedom method: satterthwaite
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 2 estimates
## P value adjustment: bonferroni method for 2 tests
summary(contra_all_P1_spec[3:4], delta = delta_null, level = .9)
summary(contra_all_P1_spec[3:4], delta = delta_null, side = ">")
summary(contra_all_P1_spec[3:4], delta = delta_null, side = "<")
Other results
Comparisons of P1 amplitudes between Durations (without multiple comparison corrections):
summary(contra_all_P1[1:8], adjust = "none")
summary(contra_all_P1[1:8], delta = delta_null, side = "two.sided", adjust = "none")
Comparisons of face-specific P1 amplitude components between Durations (without multiple comparison corrections):
summary(contra_all_P1_spec[1:4], adjust = "none")
summary(contra_all_P1_spec[1:4], delta = delta_null, side = "two.sided", adjust = "none")
Ampitudes of the N170
# only keep the data for amplitudes of the N170 (all erps)
df.erp.N170.E2 <- {
df.erp.E2.all %>%
filter(Component == "N170")
}
if (saveCSV) {
output.erp.N170.E2 <- file.path(folder_nesi, "E205_erp_N170.RData")
save(df.erp.N170.E2, file = output.erp.N170.E2)
}
# the structure of this dataset
head(df.erp.N170.E2, 10)
The maximal model
file.max.N170.E2 <- file.path(folder_nesi, "E205_N170_lmm_max.RData")
# fit the maximal model
if (!file.exists(file.max.N170.E2)) {
lmm.max.N170.E2 <- lmer(MeanAmp ~ Type * Category * Duration * Hemisphere +
(1 + Type_C + Cate_C + Dura_C + Hemi_C +
Type_Cate + Type_Dura + Cate_Dura + Type_Hemi + Cate_Hemi + Dura_Hemi +
Type_Cate_Dura + Type_Cate_Hemi + Type_Dura_Hemi + Cate_Dura_Hemi +
Type_Cate_Dura_Hemi
| SubjCode) +
(1 + Type_C + Dura_C + Hemi_C +
Type_Dura + Type_Hemi + Dura_Hemi +
Type_Dura_Hemi
| Stimuli),
data = df.erp.N170.E2,
verbose = TRUE,
control = lmerControl(optimizer = "bobyqa", # nloptwrap Nelder_Mead
optCtrl = list(maxfun = 1e7)))
} else {
load(file.max.N170.E2)
}
print(summary(lmm.max.N170.E2), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Type * Category * Duration * Hemisphere + (1 + Type_C + Cate_C + Dura_C + Hemi_C + Type_Cate + Type_Dura + Cate_Dura + Type_Hemi + Cate_Hemi + Dura_Hemi + Type_Cate_Dura + Type_Cate_Hemi + Type_Dura_Hemi + Cate_Dura_Hemi + Type_Cate_Dura_Hemi | SubjCode) + (1 + Type_C + Dura_C + Hemi_C + Type_Dura + Type_Hemi + Dura_Hemi + Type_Dura_Hemi | Stimuli)
## Data: df.erp.N170.E2
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+07))
##
## REML criterion at convergence: 401341.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.9374 -0.5879 -0.0046 0.5793 10.5966
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Stimuli (Intercept) 0.017355 0.1317
## Type_C 0.099516 0.3155 0.23
## Dura_C 0.031794 0.1783 0.62 -0.33
## Hemi_C 0.004382 0.0662 0.32 0.31 -0.14
## Type_Dura 0.375660 0.6129 0.30 0.42 0.13 -0.61
## Type_Hemi 0.002172 0.0466 0.53 0.32 0.11 -0.40 0.89
## Dura_Hemi 0.010231 0.1011 0.49 0.63 0.04 0.87 -0.25 -0.19
## Type_Dura_Hemi 0.040130 0.2003 0.17 -0.51 0.38 -0.76 0.55 0.63 -0.77
## SubjCode (Intercept) 3.125355 1.7679
## Type_C 0.509356 0.7137 -0.23
## Cate_C 0.352785 0.5940 -0.27 0.64
## Dura_C 0.906309 0.9520 0.00 0.20 0.02
## Hemi_C 1.175660 1.0843 -0.21 -0.35 -0.02 -0.34
## Type_Cate 1.443237 1.2013 0.62 -0.57 -0.79 -0.01 -0.06
## Type_Dura 0.838114 0.9155 -0.04 0.59 0.36 -0.03 -0.53 -0.21
## Cate_Dura 0.441519 0.6645 0.21 0.22 0.53 -0.20 -0.16 -0.35 0.62
## Type_Hemi 0.520714 0.7216 -0.03 0.52 0.27 0.24 -0.59 -0.27 0.43 0.21
## Cate_Hemi 0.471528 0.6867 -0.02 0.07 -0.06 -0.01 -0.42 0.02 0.24 0.15 0.83
## Dura_Hemi 1.453583 1.2056 -0.15 0.03 -0.08 0.49 -0.18 0.01 0.17 0.01 0.31 0.26
## Type_Cate_Dura 1.523837 1.2344 0.39 -0.57 -0.70 -0.14 0.22 0.75 -0.62 -0.70 -0.49 -0.30 -0.14
## Type_Cate_Hemi 1.545688 1.2433 0.04 -0.20 0.00 -0.11 0.42 0.08 -0.10 -0.03 -0.86 -0.95 -0.19 0.30
## Type_Dura_Hemi 0.518336 0.7200 0.11 0.71 0.37 -0.08 -0.48 -0.15 0.60 0.21 0.65 0.41 -0.30 -0.34 -0.45
## Cate_Dura_Hemi 0.442789 0.6654 0.11 -0.17 -0.24 -0.43 -0.08 0.21 0.22 0.24 0.48 0.79 -0.08 -0.13 -0.66 0.37
## Type_Cate_Dura_Hemi 0.644599 0.8029 0.10 -0.36 0.05 0.29 0.24 0.00 -0.46 -0.11 -0.63 -0.70 -0.02 0.23 0.65 -0.66 -0.78
## Residual 16.105886 4.0132
## Number of obs: 71282, groups: Stimuli, 80; SubjCode, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.457986 0.323589 29.121560 -4.506 9.92e-05 ***
## Type2-1 1.035785 0.139571 32.818293 7.421 1.64e-08 ***
## Category2-1 0.953981 0.117833 32.584341 8.096 2.65e-09 ***
## Duration2-1 0.560394 0.178509 29.688145 3.139 0.00381 **
## Hemisphere2-1 0.002074 0.201247 29.015948 0.010 0.99185
## Type2-1:Category2-1 -1.932099 0.241068 34.283041 -8.015 2.29e-09 ***
## Type2-1:Duration2-1 1.014670 0.194070 36.684930 5.228 7.11e-06 ***
## Category2-1:Duration2-1 1.064759 0.146072 32.945035 7.289 2.32e-08 ***
## Type2-1:Hemisphere2-1 0.388565 0.149700 29.423292 2.596 0.01458 *
## Category2-1:Hemisphere2-1 0.188786 0.144784 30.757300 1.304 0.20194
## Duration2-1:Hemisphere2-1 -0.234009 0.231533 29.172466 -1.011 0.32047
## Type2-1:Category2-1:Duration2-1 -2.037947 0.299494 42.540100 -6.805 2.60e-08 ***
## Type2-1:Category2-1:Hemisphere2-1 -0.586785 0.267832 30.765442 -2.191 0.03616 *
## Type2-1:Duration2-1:Hemisphere2-1 0.454275 0.194635 38.826366 2.334 0.02487 *
## Category2-1:Duration2-1:Hemisphere2-1 0.046010 0.188081 41.085375 0.245 0.80796
## Type2-1:Category2-1:Duration2-1:Hemisphere2-1 -0.529665 0.322344 53.344936 -1.643 0.10623
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## optimizer (bobyqa) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
The zero-correlation-parameter model
file.zcp.N170.E2 <- file.path(folder_nesi, "E205_N170_lmm_zcp.RData")
# fit the zcp model
if (!file.exists(file.zcp.N170.E2)) {
lmm.zcp.N170.E2 <- lmer(MeanAmp ~ Type * Category * Duration * Hemisphere +
(1 + Type_C + Cate_C + Dura_C + Hemi_C +
Type_Cate + Type_Dura + Cate_Dura + Type_Hemi + Cate_Hemi + Dura_Hemi +
Type_Cate_Dura + Type_Cate_Hemi + Type_Dura_Hemi + Cate_Dura_Hemi +
Type_Cate_Dura_Hemi
|| SubjCode) +
(1 + Type_C + Dura_C + Hemi_C +
Type_Dura + Type_Hemi + Dura_Hemi +
Type_Dura_Hemi
|| Stimuli),
data = df.erp.N170.E2,
verbose = TRUE,
control = lmerControl(optimizer = "bobyqa", # nloptwrap Nelder_Mead
optCtrl = list(maxfun = 1e7)))
} else {
load(file.zcp.N170.E2)
}
print(summary(lmm.zcp.N170.E2), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Type * Category * Duration * Hemisphere + (1 + Type_C + Cate_C + Dura_C + Hemi_C + Type_Cate + Type_Dura + Cate_Dura + Type_Hemi + Cate_Hemi + Dura_Hemi + Type_Cate_Dura + Type_Cate_Hemi + Type_Dura_Hemi + Cate_Dura_Hemi + Type_Cate_Dura_Hemi || SubjCode) + (1 + Type_C + Dura_C + Hemi_C + Type_Dura + Type_Hemi + Dura_Hemi + Type_Dura_Hemi || Stimuli)
## Data: df.erp.N170.E2
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+07))
##
## REML criterion at convergence: 401644.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.9367 -0.5852 -0.0038 0.5801 10.5465
##
## Random effects:
## Groups Name Variance Std.Dev.
## Stimuli Type_Dura_Hemi 1.275e-13 3.571e-07
## Stimuli.1 Dura_Hemi 2.294e-14 1.515e-07
## Stimuli.2 Type_Hemi 0.000e+00 0.000e+00
## Stimuli.3 Type_Dura 2.851e-01 5.340e-01
## Stimuli.4 Hemi_C 0.000e+00 0.000e+00
## Stimuli.5 Dura_C 8.721e-03 9.339e-02
## Stimuli.6 Type_C 6.560e-02 2.561e-01
## Stimuli.7 (Intercept) 1.123e-02 1.060e-01
## SubjCode Type_Cate_Dura_Hemi 1.072e-13 3.274e-07
## SubjCode.1 Cate_Dura_Hemi 3.997e-03 6.322e-02
## SubjCode.2 Type_Dura_Hemi 1.168e-01 3.417e-01
## SubjCode.3 Type_Cate_Hemi 1.076e+00 1.037e+00
## SubjCode.4 Type_Cate_Dura 1.287e+00 1.134e+00
## SubjCode.5 Dura_Hemi 1.510e+00 1.229e+00
## SubjCode.6 Cate_Hemi 2.930e-01 5.413e-01
## SubjCode.7 Type_Hemi 3.862e-01 6.215e-01
## SubjCode.8 Cate_Dura 4.000e-01 6.325e-01
## SubjCode.9 Type_Dura 7.758e-01 8.808e-01
## SubjCode.10 Type_Cate 1.289e+00 1.135e+00
## SubjCode.11 Hemi_C 1.186e+00 1.089e+00
## SubjCode.12 Dura_C 9.083e-01 9.530e-01
## SubjCode.13 Cate_C 3.241e-01 5.693e-01
## SubjCode.14 Type_C 4.811e-01 6.936e-01
## SubjCode.15 (Intercept) 3.128e+00 1.769e+00
## Residual 1.614e+01 4.017e+00
## Number of obs: 71282, groups: Stimuli, 80; SubjCode, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.458e+00 3.236e-01 2.908e+01 -4.506 9.93e-05 ***
## Type2-1 1.034e+00 1.346e-01 3.173e+01 7.682 9.88e-09 ***
## Category2-1 9.531e-01 1.124e-01 3.175e+01 8.483 1.15e-09 ***
## Duration2-1 5.603e-01 1.779e-01 2.920e+01 3.150 0.00376 **
## Hemisphere2-1 1.826e-03 2.020e-01 2.904e+01 0.009 0.99285
## Type2-1:Category2-1 -1.930e+00 2.264e-01 3.295e+01 -8.524 7.62e-10 ***
## Type2-1:Duration2-1 1.012e+00 1.856e-01 3.524e+01 5.450 4.02e-06 ***
## Category2-1:Duration2-1 1.063e+00 1.371e-01 2.985e+01 7.748 1.25e-08 ***
## Type2-1:Hemisphere2-1 3.885e-01 1.338e-01 3.142e+01 2.903 0.00671 **
## Category2-1:Hemisphere2-1 1.895e-01 1.217e-01 3.280e+01 1.558 0.12891
## Duration2-1:Hemisphere2-1 -2.345e-01 2.353e-01 2.914e+01 -0.996 0.32721
## Type2-1:Category2-1:Duration2-1 -2.036e+00 2.780e-01 4.017e+01 -7.322 6.54e-09 ***
## Type2-1:Category2-1:Hemisphere2-1 -5.890e-01 2.366e-01 3.446e+01 -2.489 0.01780 *
## Type2-1:Duration2-1:Hemisphere2-1 4.534e-01 1.550e-01 3.135e+01 2.925 0.00636 **
## Category2-1:Duration2-1:Hemisphere2-1 4.961e-02 1.424e-01 3.306e+01 0.348 0.72976
## Type2-1:Category2-1:Duration2-1:Hemisphere2-1 -5.319e-01 2.838e-01 7.052e+04 -1.874 0.06095 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## optimizer (bobyqa) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
The reduced model
summary(rePCA(lmm.zcp.N170.E2))
## $Stimuli
## Importance of components:
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## Standard deviation 0.1329 0.06376 0.02639 0.02325 8.89e-08 3.77e-08 0 0
## Proportion of Variance 0.7692 0.17698 0.03031 0.02353 0.00e+00 0.00e+00 0 0
## Cumulative Proportion 0.7692 0.94617 0.97647 1.00000 1.00e+00 1.00e+00 1 1
##
## $SubjCode
## Importance of components:
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16]
## Standard deviation 0.4403 0.3059 0.28258 0.28236 0.27113 0.25820 0.23725 0.21927 0.17267 0.15744 0.15471 0.14172 0.13475 0.08507 0.01574 8.15e-08
## Proportion of Variance 0.2376 0.1147 0.09788 0.09773 0.09011 0.08172 0.06899 0.05893 0.03655 0.03038 0.02934 0.02462 0.02226 0.00887 0.00030 0.00e+00
## Cumulative Proportion 0.2376 0.3523 0.45020 0.54793 0.63804 0.71975 0.78875 0.84768 0.88423 0.91461 0.94395 0.96857 0.99083 0.99970 1.00000 1.00e+00
file.rdc.N170.E2 <- file.path(folder_nesi, "E205_N170_lmm_rdc.RData")
# fit the rdc model
if (!file.exists(file.rdc.N170.E2)) {
# lmm.rdc.N170.E2.step <- step(lmm.rdc.N170.E2, reduce.fixed = FALSE)
lmm.rdc.N170.E2 <- update(lmm.zcp.N170.E2,
formula = MeanAmp ~ Type * Category * Duration * Hemisphere +
(1 + Type_C + Cate_C + Dura_C + Hemi_C +
Type_Cate + Type_Dura + Cate_Dura + Type_Hemi + Cate_Hemi + Dura_Hemi +
Type_Cate_Dura + Type_Cate_Hemi + Type_Dura_Hemi
|| SubjCode) +
(1 + Type_C + Dura_C +
Type_Dura
|| Stimuli))
} else {
load(file.rdc.N170.E2)
}
print(summary(lmm.rdc.N170.E2), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Type + Category + Duration + Hemisphere + (1 + Type_C + Cate_C + Dura_C + Hemi_C + Type_Cate + Type_Dura + Cate_Dura + Type_Hemi + Cate_Hemi + Dura_Hemi + Type_Cate_Dura + Type_Cate_Hemi + Type_Dura_Hemi || SubjCode) + (1 + Type_C + Dura_C + Type_Dura || Stimuli) + Type:Category + Type:Duration + Category:Duration + Type:Hemisphere + Category:Hemisphere + Duration:Hemisphere + Type:Category:Duration + Type:Category:Hemisphere + Type:Duration:Hemisphere + Category:Duration:Hemisphere + Type:Category:Duration:Hemisphere
## Data: df.erp.N170.E2
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+07))
##
## REML criterion at convergence: 401644.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.9367 -0.5852 -0.0037 0.5800 10.5465
##
## Random effects:
## Groups Name Variance Std.Dev.
## Stimuli Type_Dura 0.285125 0.53397
## Stimuli.1 Dura_C 0.008721 0.09339
## Stimuli.2 Type_C 0.065603 0.25613
## Stimuli.3 (Intercept) 0.011234 0.10599
## SubjCode Type_Dura_Hemi 0.116720 0.34164
## SubjCode.1 Type_Cate_Hemi 1.075684 1.03715
## SubjCode.2 Type_Cate_Dura 1.286463 1.13422
## SubjCode.3 Dura_Hemi 1.510497 1.22902
## SubjCode.4 Cate_Hemi 0.292453 0.54079
## SubjCode.5 Type_Hemi 0.386211 0.62146
## SubjCode.6 Cate_Dura 0.399993 0.63245
## SubjCode.7 Type_Dura 0.775818 0.88081
## SubjCode.8 Type_Cate 1.288519 1.13513
## SubjCode.9 Hemi_C 1.186257 1.08915
## SubjCode.10 Dura_C 0.908297 0.95305
## SubjCode.11 Cate_C 0.324088 0.56929
## SubjCode.12 Type_C 0.481105 0.69362
## SubjCode.13 (Intercept) 3.127930 1.76860
## Residual 16.137008 4.01709
## Number of obs: 71282, groups: Stimuli, 80; SubjCode, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.458e+00 3.236e-01 2.908e+01 -4.506 9.93e-05 ***
## Type2-1 1.034e+00 1.346e-01 3.173e+01 7.682 9.88e-09 ***
## Category2-1 9.531e-01 1.124e-01 3.175e+01 8.483 1.15e-09 ***
## Duration2-1 5.603e-01 1.779e-01 2.920e+01 3.150 0.00376 **
## Hemisphere2-1 1.825e-03 2.020e-01 2.904e+01 0.009 0.99285
## Type2-1:Category2-1 -1.930e+00 2.264e-01 3.294e+01 -8.524 7.62e-10 ***
## Type2-1:Duration2-1 1.012e+00 1.856e-01 3.524e+01 5.450 4.02e-06 ***
## Category2-1:Duration2-1 1.063e+00 1.371e-01 2.985e+01 7.748 1.25e-08 ***
## Type2-1:Hemisphere2-1 3.885e-01 1.338e-01 3.143e+01 2.903 0.00671 **
## Category2-1:Hemisphere2-1 1.895e-01 1.216e-01 3.413e+01 1.559 0.12829
## Duration2-1:Hemisphere2-1 -2.345e-01 2.353e-01 2.914e+01 -0.996 0.32722
## Type2-1:Category2-1:Duration2-1 -2.036e+00 2.780e-01 4.017e+01 -7.322 6.54e-09 ***
## Type2-1:Category2-1:Hemisphere2-1 -5.890e-01 2.366e-01 3.446e+01 -2.489 0.01780 *
## Type2-1:Duration2-1:Hemisphere2-1 4.534e-01 1.550e-01 3.135e+01 2.925 0.00636 **
## Category2-1:Duration2-1:Hemisphere2-1 4.963e-02 1.419e-01 7.054e+04 0.350 0.72655
## Type2-1:Category2-1:Duration2-1:Hemisphere2-1 -5.319e-01 2.838e-01 7.054e+04 -1.874 0.06093 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lmm.rdc.N170.E2, lmm.zcp.N170.E2, refit = FALSE)
The extended model
file.etd.N170.E2 <- file.path(folder_nesi, "E205_N170_lmm_etd.RData")
# fit the etd model
if (!file.exists(file.etd.N170.E2)) {
lmm.etd.N170.E2 <- update(lmm.rdc.N170.E2,
formula = MeanAmp ~ Type * Category * Duration * Hemisphere +
(1 + Type_C + Cate_C + Dura_C + Hemi_C +
Type_Cate + Type_Dura + Cate_Dura + Type_Hemi + Cate_Hemi + Dura_Hemi +
Type_Cate_Dura + Type_Cate_Hemi + Type_Dura_Hemi
| SubjCode) +
(1 + Type_C + Dura_C +
Type_Dura
| Stimuli))
} else {
load(file.etd.N170.E2)
}
print(summary(lmm.etd.N170.E2), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Type + Category + Duration + Hemisphere + (1 + Type_C + Cate_C + Dura_C + Hemi_C + Type_Cate + Type_Dura + Cate_Dura + Type_Hemi + Cate_Hemi + Dura_Hemi + Type_Cate_Dura + Type_Cate_Hemi + Type_Dura_Hemi | SubjCode) + (1 + Type_C + Dura_C + Type_Dura | Stimuli) + Type:Category + Type:Duration + Category:Duration + Type:Hemisphere + Category:Hemisphere + Duration:Hemisphere + Type:Category:Duration + Type:Category:Hemisphere + Type:Duration:Hemisphere + Category:Duration:Hemisphere + Type:Category:Duration:Hemisphere
## Data: df.erp.N170.E2
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+07))
##
## REML criterion at convergence: 401372.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.9468 -0.5869 -0.0042 0.5788 10.5930
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Stimuli (Intercept) 0.01700 0.1304
## Type_C 0.09751 0.3123 0.22
## Dura_C 0.03059 0.1749 0.63 -0.36
## Type_Dura 0.36712 0.6059 0.32 0.43 0.12
## SubjCode (Intercept) 3.12552 1.7679
## Type_C 0.50943 0.7137 -0.23
## Cate_C 0.35290 0.5941 -0.27 0.64
## Dura_C 0.90612 0.9519 0.00 0.20 0.02
## Hemi_C 1.17515 1.0840 -0.20 -0.35 -0.02 -0.34
## Type_Cate 1.44334 1.2014 0.62 -0.57 -0.79 -0.01 -0.06
## Type_Dura 0.83973 0.9164 -0.04 0.59 0.36 -0.03 -0.53 -0.21
## Cate_Dura 0.44172 0.6646 0.21 0.21 0.53 -0.20 -0.16 -0.35 0.61
## Type_Hemi 0.51510 0.7177 -0.03 0.53 0.27 0.24 -0.60 -0.28 0.44 0.20
## Cate_Hemi 0.30258 0.5501 -0.05 0.12 0.00 0.12 -0.49 -0.04 0.21 0.11 0.86
## Dura_Hemi 1.45408 1.2059 -0.15 0.03 -0.08 0.49 -0.18 0.01 0.17 0.01 0.31 0.35
## Type_Cate_Dura 1.52322 1.2342 0.39 -0.58 -0.70 -0.14 0.22 0.75 -0.62 -0.70 -0.50 -0.33 -0.14
## Type_Cate_Hemi 1.16352 1.0787 0.04 -0.18 -0.02 -0.22 0.43 0.11 -0.02 0.01 -0.85 -0.95 -0.22 0.29
## Type_Dura_Hemi 0.49472 0.7034 0.11 0.72 0.39 -0.08 -0.49 -0.16 0.61 0.20 0.64 0.37 -0.31 -0.35 -0.37
## Residual 16.11624 4.0145
## Number of obs: 71282, groups: Stimuli, 80; SubjCode, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.458e+00 3.236e-01 2.900e+01 -4.506 9.99e-05 ***
## Type2-1 1.036e+00 1.395e-01 2.389e+01 7.425 1.19e-07 ***
## Category2-1 9.539e-01 1.178e-01 2.705e+01 8.099 1.05e-08 ***
## Duration2-1 5.603e-01 1.785e-01 2.390e+01 3.140 0.00445 **
## Hemisphere2-1 1.499e-03 2.011e-01 1.289e+01 0.007 0.99417
## Type2-1:Category2-1 -1.932e+00 2.409e-01 2.846e+01 -8.021 8.72e-09 ***
## Type2-1:Duration2-1 1.015e+00 1.939e-01 3.399e+01 5.232 8.60e-06 ***
## Category2-1:Duration2-1 1.065e+00 1.459e-01 3.290e+01 7.298 2.29e-08 ***
## Type2-1:Hemisphere2-1 3.889e-01 1.490e-01 2.525e+01 2.610 0.01501 *
## Category2-1:Hemisphere2-1 1.901e-01 1.229e-01 3.634e+01 1.546 0.13075
## Duration2-1:Hemisphere2-1 -2.345e-01 2.313e-01 1.980e+01 -1.014 0.32287
## Type2-1:Category2-1:Duration2-1 -2.038e+00 2.988e-01 3.953e+01 -6.821 3.53e-08 ***
## Type2-1:Category2-1:Hemisphere2-1 -5.894e-01 2.427e-01 3.582e+01 -2.429 0.02031 *
## Type2-1:Duration2-1:Hemisphere2-1 4.561e-01 1.913e-01 3.071e+01 2.384 0.02351 *
## Category2-1:Duration2-1:Hemisphere2-1 4.947e-02 1.418e-01 7.067e+04 0.349 0.72720
## Type2-1:Category2-1:Duration2-1:Hemisphere2-1 -5.344e-01 2.836e-01 7.067e+04 -1.884 0.05956 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## optimizer (bobyqa) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
summary(rePCA(lmm.etd.N170.E2))
## $Stimuli
## Importance of components:
## [,1] [,2] [,3] [,4]
## Standard deviation 0.1559 0.07239 0.04451 0.01604
## Proportion of Variance 0.7647 0.16489 0.06235 0.00809
## Cumulative Proportion 0.7647 0.92956 0.99191 1.00000
##
## $SubjCode
## Importance of components:
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
## Standard deviation 0.5573 0.4522 0.3638 0.28544 0.23566 0.21812 0.16008 0.13781 0.06854 0.05289 0.0001626 2.121e-06 6.643e-08 0
## Proportion of Variance 0.3513 0.2313 0.1497 0.09216 0.06282 0.05382 0.02899 0.02148 0.00531 0.00316 0.0000000 0.000e+00 0.000e+00 0
## Cumulative Proportion 0.3513 0.5826 0.7322 0.82441 0.88724 0.94105 0.97004 0.99152 0.99684 1.00000 1.0000000 1.000e+00 1.000e+00 1
file.etd1.N170.E2 <- file.path(folder_nesi, "E205_N170_lmm_etd1.RData")
# fit the etd1 model
if (!file.exists(file.etd1.N170.E2)) {
lmm.etd1.N170.E2 <- update(lmm.etd.N170.E2,
formula = MeanAmp ~ Type * Category * Duration * Hemisphere +
(1 + Type_C + Dura_C + Hemi_C +
Type_Cate + Type_Dura + Type_Hemi + Dura_Hemi +
Type_Cate_Dura + Type_Cate_Hemi
| SubjCode) +
(1 + Type_C + Dura_C +
Type_Dura
| Stimuli))
} else {
load(file.etd1.N170.E2)
}
print(summary(lmm.etd1.N170.E2), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Type + Category + Duration + Hemisphere + (1 + Type_C + Dura_C + Hemi_C + Type_Cate + Type_Dura + Type_Hemi + Dura_Hemi + Type_Cate_Dura + Type_Cate_Hemi | SubjCode) + (1 + Type_C + Dura_C + Type_Dura | Stimuli) + Type:Category + Type:Duration + Category:Duration + Type:Hemisphere + Category:Hemisphere + Duration:Hemisphere + Type:Category:Duration + Type:Category:Hemisphere + Type:Duration:Hemisphere + Category:Duration:Hemisphere + Type:Category:Duration:Hemisphere
## Data: df.erp.N170.E2
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+07))
##
## REML criterion at convergence: 401774.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.9501 -0.5858 -0.0039 0.5816 10.6037
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Stimuli (Intercept) 0.01651 0.1285
## Type_C 0.09623 0.3102 0.21
## Dura_C 0.03051 0.1747 0.64 -0.36
## Type_Dura 0.36483 0.6040 0.33 0.42 0.13
## SubjCode (Intercept) 3.12440 1.7676
## Type_C 0.50975 0.7140 -0.23
## Dura_C 0.90542 0.9515 0.00 0.20
## Hemi_C 1.15652 1.0754 -0.21 -0.34 -0.34
## Type_Cate 1.75252 1.3238 0.62 -0.59 -0.03 -0.07
## Type_Dura 0.84363 0.9185 -0.04 0.60 -0.03 -0.53 -0.20
## Type_Hemi 0.38161 0.6177 -0.06 0.42 0.30 -0.56 -0.28 0.34
## Dura_Hemi 1.48018 1.2166 -0.15 0.01 0.49 -0.17 0.03 0.15 0.43
## Type_Cate_Dura 1.31094 1.1450 0.30 -0.44 -0.10 0.25 0.55 -0.66 -0.49 -0.17
## Type_Cate_Hemi 1.69344 1.3013 0.04 -0.18 -0.21 0.46 0.09 -0.05 -0.90 -0.25 0.33
## Residual 16.22266 4.0277
## Number of obs: 71282, groups: Stimuli, 80; SubjCode, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.457e+00 3.235e-01 2.912e+01 -4.504 9.97e-05 ***
## Type2-1 1.035e+00 1.395e-01 3.271e+01 7.417 1.69e-08 ***
## Category2-1 9.523e-01 4.574e-02 7.773e+01 20.822 < 2e-16 ***
## Duration2-1 5.637e-01 1.784e-01 2.968e+01 3.160 0.00362 **
## Hemisphere2-1 2.036e-03 1.995e-01 2.905e+01 0.010 0.99193
## Type2-1:Category2-1 -1.928e+00 2.613e-01 3.344e+01 -7.379 1.64e-08 ***
## Type2-1:Duration2-1 1.017e+00 1.943e-01 3.637e+01 5.235 7.13e-06 ***
## Category2-1:Duration2-1 1.060e+00 8.118e-02 7.754e+01 13.057 < 2e-16 ***
## Type2-1:Hemisphere2-1 3.891e-01 1.333e-01 3.398e+01 2.918 0.00620 **
## Category2-1:Hemisphere2-1 1.909e-01 7.114e-02 7.070e+04 2.683 0.00730 **
## Duration2-1:Hemisphere2-1 -2.350e-01 2.332e-01 2.916e+01 -1.007 0.32199
## Type2-1:Category2-1:Duration2-1 -2.031e+00 2.867e-01 4.327e+01 -7.086 9.37e-09 ***
## Type2-1:Category2-1:Hemisphere2-1 -5.876e-01 2.769e-01 3.445e+01 -2.122 0.04111 *
## Type2-1:Duration2-1:Hemisphere2-1 4.548e-01 1.423e-01 7.070e+04 3.196 0.00139 **
## Category2-1:Duration2-1:Hemisphere2-1 4.671e-02 1.423e-01 7.070e+04 0.328 0.74268
## Type2-1:Category2-1:Duration2-1:Hemisphere2-1 -5.362e-01 2.846e-01 7.070e+04 -1.884 0.05956 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## optimizer (bobyqa) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
summary(rePCA(lmm.etd1.N170.E2))
## $Stimuli
## Importance of components:
## [,1] [,2] [,3] [,4]
## Standard deviation 0.1548 0.07198 0.04396 0.01573
## Proportion of Variance 0.7650 0.16543 0.06171 0.00790
## Cumulative Proportion 0.7650 0.93039 0.99210 1.00000
##
## $SubjCode
## Importance of components:
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## Standard deviation 0.5368 0.4428 0.3314 0.28204 0.22279 0.21459 0.1596 0.12076 0.04202 6.325e-08
## Proportion of Variance 0.3553 0.2418 0.1354 0.09807 0.06119 0.05677 0.0314 0.01798 0.00218 0.000e+00
## Cumulative Proportion 0.3553 0.5970 0.7324 0.83048 0.89168 0.94845 0.9798 0.99782 1.00000 1.000e+00
file.etd2.N170.E2 <- file.path(folder_nesi, "E205_N170_lmm_etd2.RData")
# fit the etd2 model
if (!file.exists(file.etd2.N170.E2)) {
lmm.etd2.N170.E2 <- update(lmm.etd1.N170.E2,
formula = MeanAmp ~ Type * Category * Duration * Hemisphere +
(1 + Type_C + Dura_C + Hemi_C +
Type_Cate + Type_Dura + Dura_Hemi +
Type_Cate_Dura + Type_Cate_Hemi
| SubjCode) +
(1 + Type_C + Dura_C +
Type_Dura
| Stimuli))
} else {
load(file.etd2.N170.E2)
}
print(summary(lmm.etd2.N170.E2), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Type + Category + Duration + Hemisphere + (1 + Type_C + Dura_C + Hemi_C + Type_Cate + Type_Dura + Dura_Hemi + Type_Cate_Dura + Type_Cate_Hemi | SubjCode) + (1 + Type_C + Dura_C + Type_Dura | Stimuli) + Type:Category + Type:Duration + Category:Duration + Type:Hemisphere + Category:Hemisphere + Duration:Hemisphere + Type:Category:Duration + Type:Category:Hemisphere + Type:Duration:Hemisphere + Category:Duration:Hemisphere + Type:Category:Duration:Hemisphere
## Data: df.erp.N170.E2
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+07))
##
## REML criterion at convergence: 401871.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.9303 -0.5860 -0.0048 0.5817 10.5986
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Stimuli (Intercept) 0.002805 0.05297
## Type_C 0.092295 0.30380 1.00
## Dura_C 0.011615 0.10777 -0.32 -0.32
## Type_Dura 0.399586 0.63213 0.44 0.44 0.08
## SubjCode (Intercept) 3.123796 1.76743
## Type_C 0.509423 0.71374 -0.23
## Dura_C 0.905459 0.95156 0.00 0.20
## Hemi_C 1.218893 1.10403 -0.20 -0.35 -0.35
## Type_Cate 1.753881 1.32434 0.62 -0.59 -0.03 -0.05
## Type_Dura 0.838924 0.91593 -0.04 0.60 -0.03 -0.53 -0.20
## Dura_Hemi 1.601393 1.26546 -0.15 0.04 0.49 -0.23 0.00 0.17
## Type_Cate_Dura 1.313151 1.14593 0.30 -0.44 -0.10 0.27 0.55 -0.67 -0.21
## Type_Cate_Hemi 1.594705 1.26282 0.04 -0.18 -0.21 0.50 0.09 -0.04 -0.31 0.36
## Residual 16.256515 4.03194
## Number of obs: 71282, groups: Stimuli, 80; SubjCode, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.457e+00 3.232e-01 2.902e+01 -4.507 9.93e-05 ***
## Type2-1 1.034e+00 1.393e-01 3.248e+01 7.425 1.73e-08 ***
## Category2-1 9.525e-01 3.754e-02 4.821e+01 25.376 < 2e-16 ***
## Duration2-1 5.641e-01 1.778e-01 2.926e+01 3.174 0.00353 **
## Hemisphere2-1 2.941e-03 2.047e-01 2.905e+01 0.014 0.98864
## Type2-1:Category2-1 -1.929e+00 2.611e-01 3.316e+01 -7.387 1.69e-08 ***
## Type2-1:Duration2-1 1.016e+00 1.950e-01 3.715e+01 5.208 7.32e-06 ***
## Category2-1:Duration2-1 1.060e+00 7.520e-02 1.147e+02 14.100 < 2e-16 ***
## Type2-1:Hemisphere2-1 3.885e-01 7.122e-02 7.080e+04 5.455 4.92e-08 ***
## Category2-1:Hemisphere2-1 1.916e-01 7.122e-02 7.080e+04 2.690 0.00714 **
## Duration2-1:Hemisphere2-1 -2.349e-01 2.418e-01 2.922e+01 -0.972 0.33917
## Type2-1:Category2-1:Duration2-1 -2.031e+00 2.899e-01 4.504e+01 -7.006 9.93e-09 ***
## Type2-1:Category2-1:Hemisphere2-1 -5.857e-01 2.710e-01 3.421e+01 -2.161 0.03776 *
## Type2-1:Duration2-1:Hemisphere2-1 4.493e-01 1.424e-01 7.080e+04 3.155 0.00161 **
## Category2-1:Duration2-1:Hemisphere2-1 4.771e-02 1.424e-01 7.080e+04 0.335 0.73767
## Type2-1:Category2-1:Duration2-1:Hemisphere2-1 -5.345e-01 2.849e-01 7.080e+04 -1.876 0.06065 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## optimizer (bobyqa) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
summary(rePCA(lmm.etd2.N170.E2))
## $Stimuli
## Importance of components:
## [,1] [,2] [,3] [,4]
## Standard deviation 0.1611 0.0678 0.02407 1.847e-18
## Proportion of Variance 0.8338 0.1476 0.01861 0.000e+00
## Cumulative Proportion 0.8338 0.9814 1.00000 1.000e+00
##
## $SubjCode
## Importance of components:
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## Standard deviation 0.5344 0.4352 0.3311 0.26908 0.2225 0.21374 0.15795 0.11781 2.884e-05
## Proportion of Variance 0.3610 0.2394 0.1386 0.09153 0.0626 0.05775 0.03154 0.01755 0.000e+00
## Cumulative Proportion 0.3610 0.6004 0.7390 0.83057 0.8932 0.95092 0.98245 1.00000 1.000e+00
file.etd3.N170.E2 <- file.path(folder_nesi, "E205_N170_lmm_etd3.RData")
# fit the etd3 model
if (!file.exists(file.etd3.N170.E2)) {
lmm.etd3.N170.E2 <- update(lmm.etd2.N170.E2,
formula = MeanAmp ~ Type * Category * Duration * Hemisphere +
(1 + Dura_C + Hemi_C +
Type_Cate + Type_Dura + Dura_Hemi +
Type_Cate_Dura + Type_Cate_Hemi
| SubjCode) +
(0 + Type_C + Dura_C +
Type_Dura
| Stimuli))
} else {
load(file.etd3.N170.E2)
}
print(summary(lmm.etd3.N170.E2), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Type + Category + Duration + Hemisphere + (1 + Dura_C + Hemi_C + Type_Cate + Type_Dura + Dura_Hemi + Type_Cate_Dura + Type_Cate_Hemi | SubjCode) + (0 + Type_C + Dura_C + Type_Dura | Stimuli) + Type:Category + Type:Duration + Category:Duration + Type:Hemisphere + Category:Hemisphere + Duration:Hemisphere + Type:Category:Duration + Type:Category:Hemisphere + Type:Duration:Hemisphere + Category:Duration:Hemisphere + Type:Category:Duration:Hemisphere
## Data: df.erp.N170.E2
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+07))
##
## REML criterion at convergence: 402199
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.9564 -0.5894 -0.0039 0.5825 10.5129
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Stimuli Type_C 0.09270 0.3045
## Dura_C 0.01805 0.1343 -0.66
## Type_Dura 0.38023 0.6166 0.41 -0.14
## SubjCode (Intercept) 3.14852 1.7744
## Dura_C 0.92884 0.9638 -0.02
## Hemi_C 1.21852 1.1039 -0.19 -0.37
## Type_Cate 1.75483 1.3247 0.63 -0.07 -0.05
## Type_Dura 0.54862 0.7407 0.15 -0.23 -0.34 0.29
## Dura_Hemi 1.60033 1.2650 -0.15 0.49 -0.23 0.00 0.17
## Type_Cate_Dura 1.31892 1.1484 0.30 -0.13 0.27 0.55 -0.43 -0.20
## Type_Cate_Hemi 1.59093 1.2613 0.05 -0.22 0.50 0.09 0.12 -0.31 0.36
## Residual 16.34600 4.0430
## Number of obs: 71282, groups: Stimuli, 80; SubjCode, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.456e+00 3.245e-01 2.900e+01 -4.488 0.000105 ***
## Type2-1 1.031e+00 4.934e-02 7.879e+01 20.905 < 2e-16 ***
## Category2-1 9.524e-01 3.571e-02 7.090e+04 26.665 < 2e-16 ***
## Duration2-1 5.661e-01 1.802e-01 2.942e+01 3.142 0.003812 **
## Hemisphere2-1 2.973e-03 2.047e-01 2.905e+01 0.015 0.988510
## Type2-1:Category2-1 -1.924e+00 2.612e-01 3.328e+01 -7.367 1.76e-08 ***
## Type2-1:Duration2-1 1.011e+00 1.678e-01 4.324e+01 6.025 3.30e-07 ***
## Category2-1:Duration2-1 1.060e+00 7.749e-02 1.154e+02 13.685 < 2e-16 ***
## Type2-1:Hemisphere2-1 3.885e-01 7.142e-02 7.083e+04 5.439 5.36e-08 ***
## Category2-1:Hemisphere2-1 1.916e-01 7.141e-02 7.083e+04 2.683 0.007308 **
## Duration2-1:Hemisphere2-1 -2.349e-01 2.418e-01 2.921e+01 -0.972 0.339161
## Type2-1:Category2-1:Duration2-1 -2.021e+00 2.888e-01 4.417e+01 -7.000 1.12e-08 ***
## Type2-1:Category2-1:Hemisphere2-1 -5.858e-01 2.710e-01 3.424e+01 -2.162 0.037722 *
## Type2-1:Duration2-1:Hemisphere2-1 4.492e-01 1.428e-01 7.083e+04 3.145 0.001661 **
## Category2-1:Duration2-1:Hemisphere2-1 4.772e-02 1.428e-01 7.083e+04 0.334 0.738281
## Type2-1:Category2-1:Duration2-1:Hemisphere2-1 -5.343e-01 2.857e-01 7.083e+04 -1.870 0.061423 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lmm.etd3.N170.E2, lmm.rdc.N170.E2, refit = FALSE)
The optimal model
lmm.opt.N170.E205 <- lmm.rdc.N170.E2
print(summary(lmm.opt.N170.E205), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Type + Category + Duration + Hemisphere + (1 + Type_C + Cate_C + Dura_C + Hemi_C + Type_Cate + Type_Dura + Cate_Dura + Type_Hemi + Cate_Hemi + Dura_Hemi + Type_Cate_Dura + Type_Cate_Hemi + Type_Dura_Hemi || SubjCode) + (1 + Type_C + Dura_C + Type_Dura || Stimuli) + Type:Category + Type:Duration + Category:Duration + Type:Hemisphere + Category:Hemisphere + Duration:Hemisphere + Type:Category:Duration + Type:Category:Hemisphere + Type:Duration:Hemisphere + Category:Duration:Hemisphere + Type:Category:Duration:Hemisphere
## Data: df.erp.N170.E2
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+07))
##
## REML criterion at convergence: 401644.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.9367 -0.5852 -0.0037 0.5800 10.5465
##
## Random effects:
## Groups Name Variance Std.Dev.
## Stimuli Type_Dura 0.285125 0.53397
## Stimuli.1 Dura_C 0.008721 0.09339
## Stimuli.2 Type_C 0.065603 0.25613
## Stimuli.3 (Intercept) 0.011234 0.10599
## SubjCode Type_Dura_Hemi 0.116720 0.34164
## SubjCode.1 Type_Cate_Hemi 1.075684 1.03715
## SubjCode.2 Type_Cate_Dura 1.286463 1.13422
## SubjCode.3 Dura_Hemi 1.510497 1.22902
## SubjCode.4 Cate_Hemi 0.292453 0.54079
## SubjCode.5 Type_Hemi 0.386211 0.62146
## SubjCode.6 Cate_Dura 0.399993 0.63245
## SubjCode.7 Type_Dura 0.775818 0.88081
## SubjCode.8 Type_Cate 1.288519 1.13513
## SubjCode.9 Hemi_C 1.186257 1.08915
## SubjCode.10 Dura_C 0.908297 0.95305
## SubjCode.11 Cate_C 0.324088 0.56929
## SubjCode.12 Type_C 0.481105 0.69362
## SubjCode.13 (Intercept) 3.127930 1.76860
## Residual 16.137008 4.01709
## Number of obs: 71282, groups: Stimuli, 80; SubjCode, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.458e+00 3.236e-01 2.908e+01 -4.506 9.93e-05 ***
## Type2-1 1.034e+00 1.346e-01 3.173e+01 7.682 9.88e-09 ***
## Category2-1 9.531e-01 1.124e-01 3.175e+01 8.483 1.15e-09 ***
## Duration2-1 5.603e-01 1.779e-01 2.920e+01 3.150 0.00376 **
## Hemisphere2-1 1.825e-03 2.020e-01 2.904e+01 0.009 0.99285
## Type2-1:Category2-1 -1.930e+00 2.264e-01 3.294e+01 -8.524 7.62e-10 ***
## Type2-1:Duration2-1 1.012e+00 1.856e-01 3.524e+01 5.450 4.02e-06 ***
## Category2-1:Duration2-1 1.063e+00 1.371e-01 2.985e+01 7.748 1.25e-08 ***
## Type2-1:Hemisphere2-1 3.885e-01 1.338e-01 3.143e+01 2.903 0.00671 **
## Category2-1:Hemisphere2-1 1.895e-01 1.216e-01 3.413e+01 1.559 0.12829
## Duration2-1:Hemisphere2-1 -2.345e-01 2.353e-01 2.914e+01 -0.996 0.32722
## Type2-1:Category2-1:Duration2-1 -2.036e+00 2.780e-01 4.017e+01 -7.322 6.54e-09 ***
## Type2-1:Category2-1:Hemisphere2-1 -5.890e-01 2.366e-01 3.446e+01 -2.489 0.01780 *
## Type2-1:Duration2-1:Hemisphere2-1 4.534e-01 1.550e-01 3.135e+01 2.925 0.00636 **
## Category2-1:Duration2-1:Hemisphere2-1 4.963e-02 1.419e-01 7.054e+04 0.350 0.72655
## Type2-1:Category2-1:Duration2-1:Hemisphere2-1 -5.319e-01 2.838e-01 7.054e+04 -1.874 0.06093 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Diagnostic plots
# qqplot
qqplot_lmer(lmm.opt.N170.E205)
Estimated marginal means
emm.N170.E205 <- emmeans(lmm.opt.N170.E205, ~ Type + Hemisphere + Category + Duration)
emm.N170.E205 %>%
as.data.frame() %>%
arrange(Type, Hemisphere, Category, Duration)
Plots
# emmip(emm.N170.E205, Category ~ Duration | Type + Hemisphere, CIs = TRUE) +
# scale_y_continuous(trans = "reverse")
n170.all.LinePlot = {
ggplot(data = data.frame(emm.N170.E205), aes(y = emmean, x = Duration, color = Category, group = Category)) + # set the data, varialbes for x and y axises, and the variable for dividing data
geom_point(size = 2) +
geom_line(aes(linetype = Category), size = .85) + # position = "dodge", alpha = .7
scale_linetype_manual(values=c("solid", "dashed")) +
geom_errorbar(mapping = aes(ymin = lower.CL, ymax = upper.CL), linetype = 1, # set the error bar
show.legend = FALSE, width = 0.25, alpha = .5) + # , position = position_dodge(width=0.9)
facet_grid(Type ~ Hemisphere, # switch = "x",
labeller = labeller(Hemisphere = hemiLabel)) +
scale_y_reverse(limits = ylimit_n170, breaks = seq(1, -5, -1)) + # set the limit for y axis
# scale_y_continuous(expand= c(0, 0)) + # remove the space between columns and x axis
labs(x = "Duration (ms)", y = expression(paste("Amplitudes (", mu, "V)")), fill = "Category") + # set the names for main, x and y axises title = "",
# scale_x_discrete(labels = c("left", "right")) +
# geom_text(label = c("", "", "", "", "***", "", "", ""), color = "red", size = 6, nudge_y = 0.5, nudge_x = 0.5) + # add starts to the significant columns
theme_bw() +
erp_theme
}
# n170.all.LinePlot.slide <- {
# n170.all.LinePlot +
# scale_y_reverse(limits = c(0.8, -4.8), breaks = seq(0, -4, -1)) + # set the limit for y axis
# theme(legend.position="right",
# axis.line = element_line(color="black", size = 3))
# }
# ggsave('N170_all.png', n170.all.LinePlot.slide, width = 7, height = 4)
n170.all.LinePlot
n170_all_LinePlot <- {
data.frame(emm.N170.E205) %>%
mutate(Duration = fct_recode(Duration, `33`="17", `216`="200")) %>%
ggplot(aes(y = emmean, x = Duration, group = Category)) + # set the data, varialbes for x and y axises, and the variable for dividing data
geom_point(aes(fill = Duration), size=4, shape=21, stroke=0) +
scale_fill_manual(values = duration_color[c(3, 4)], guide = "none") +
geom_line(aes(linetype = Category)) + # position = "dodge", alpha = .7
scale_linetype_manual(values=c("solid", "dashed")) +
geom_errorbar(mapping = aes(ymin = lower.CL, ymax = upper.CL), linetype = 1, # set the error bar
show.legend = FALSE, width = 0.25, alpha = .5) + # , position = position_dodge(width=0.9)
facet_grid(Type ~ Hemisphere, # switch = "x",
labeller = labeller(Hemisphere = hemiLabel)) +
scale_y_reverse(limits = ylimit_n170, breaks = seq(1, -5, -1)) + # set the limit for y axis
# scale_y_continuous(expand= c(0, 0)) + # remove the space between columns and x axis
labs(x = "Duration (ms)", y = expression(paste("N170 Amplitudes (", mu, "V)")), color = "Category", linetype = "Category") + # set the names for main, x and y axises title = "",
# scale_x_discrete(labels = c("left", "right")) +
# geom_text(label = c("", "", "", "", "***", "", "", ""), color = "red", size = 6, nudge_y = 0.5, nudge_x = 0.5) + # add starts to the significant columns
theme_bw() +
erp_theme +
theme(
legend.position = c(.5, .5),
legend.title=element_text(size=10),
legend.text=element_text(size=10),
) +
guides(linetype = guide_legend(title.position = "left", nrow = 1),
color = guide_legend(title.position = "left", nrow = 1))
}
ggsave('N170_all_pub.pdf', n170_all_LinePlot, width = 6, height = 5)
Contrasts
Main results for intact stimuli
Comparisons of N170 amplitudes for intact faces with different durations:
contra_all_N170 <- contrast(emmeans(lmm.opt.N170.E205, ~ Duration | Hemisphere + Category + Type),
"pairwise", adjust = "Bonferroni", infer = TRUE)
contra_all_N170[1:2]
## contrast Hemisphere Category Type estimate SE df lower.CL upper.CL t.ratio p.value
## 17 - 200 Left face intact 0.6763417 0.2582103 100.06 0.0887499 1.263934 2.619 0.0204
## 17 - 200 Right face intact 1.2953806 0.2582103 100.06 0.7077888 1.882972 5.017 <.0001
##
## Degrees-of-freedom method: satterthwaite
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 2 estimates
## P value adjustment: bonferroni method for 2 tests
summary(contra_all_N170[1:2], delta = delta_null, level = .9)
summary(contra_all_N170[1:2], delta = delta_null, side = ">")
summary(contra_all_N170[1:2], delta = delta_null, side = "<")
Comparisons of N170 amplitude lateralization for intact faces with different durations:
contra_all_N170_hemi <- contrast(emmeans(lmm.opt.N170.E205, ~ Duration + Hemisphere | Category + Type),
interaction = "pairwise", infer = TRUE)
contra_all_N170_hemi[1]
## Duration_pairwise Hemisphere_pairwise Category Type estimate SE df lower.CL upper.CL t.ratio p.value
## 17 - 200 Left - Right face intact -0.6190389 0.2636081 45.58 -1.149785 -0.0882927 -2.348 0.0233
##
## Degrees-of-freedom method: satterthwaite
## Confidence level used: 0.95
summary(contra_all_N170_hemi[1], delta = delta_null, level = .9)
summary(contra_all_N170_hemi[1], delta = delta_null, side = ">")
summary(contra_all_N170_hemi[1], delta = delta_null, side = "<")
Comparisons of face-specific N170 amplitude components for intact stimuli with different durations:
contra_all_N170_spec <- contrast(emmeans(lmm.opt.N170.E205, ~ Category + Duration | Hemisphere + Type),
interaction = "pairwise", adjust = "Bonferroni", infer = TRUE)
contra_all_N170_spec[1:2]
## Category_pairwise Duration_pairwise Hemisphere Type estimate SE df lower.CL upper.CL t.ratio p.value
## face - house 17 - 200 Left intact 1.922618 0.2151955 104.11 1.433203 2.412034 8.934 <.0001
## face - house 17 - 200 Right intact 2.238221 0.2151955 104.11 1.748805 2.727636 10.401 <.0001
##
## Degrees-of-freedom method: satterthwaite
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 2 estimates
## P value adjustment: bonferroni method for 2 tests
summary(contra_all_N170_spec[1:2], delta = delta_null, level = .9)
summary(contra_all_N170_spec[1:2], delta = delta_null, side = ">")
summary(contra_all_N170_spec[1:2], delta = delta_null, side = "<")
Comparisons of face-specific N170 amplitude component lateralization for intact stimuli with different durations:
contra_all_N170_spec_hemi <- contrast(emmeans(lmm.opt.N170.E205, ~ Category + Duration + Hemisphere | Type),
interaction = "pairwise", adjust = "Bonferroni", infer = TRUE)
contra_all_N170_spec_hemi[1]
## Category_pairwise Duration_pairwise Hemisphere_pairwise Type estimate SE df lower.CL upper.CL t.ratio p.value
## face - house 17 - 200 Left - Right intact -0.3156022 0.1910371 70540.74 -0.6900345 0.05883006 -1.652 0.0985
##
## Degrees-of-freedom method: satterthwaite
## Confidence level used: 0.95
summary(contra_all_N170_spec_hemi[1], delta = delta_null, level = .9)
summary(contra_all_N170_spec_hemi[1], delta = delta_null, side = ">")
summary(contra_all_N170_spec_hemi[1], delta = delta_null, side = "<")
Main results for scrambled stimuli
Comparisons of N170 amplitudes for scrambled faces with different durations:
contra_all_N170[5:6]
## contrast Hemisphere Category Type estimate SE df lower.CL upper.CL t.ratio p.value
## 17 - 200 Left face scrambled -0.9936544 0.2619497 105.98 -1.589248 -0.3980611 -3.793 0.0005
## 17 - 200 Right face scrambled -1.0940310 0.2619497 105.98 -1.689624 -0.4984377 -4.176 0.0001
##
## Degrees-of-freedom method: satterthwaite
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 2 estimates
## P value adjustment: bonferroni method for 2 tests
summary(contra_all_N170[5:6], delta = delta_null, level = .9)
summary(contra_all_N170[5:6], delta = delta_null, side = ">")
summary(contra_all_N170[5:6], delta = delta_null, side = "<")
Comparisons of N170 amplitude lateralization for scrambled faces with different durations:
contra_all_N170_hemi[3]
## Duration_pairwise Hemisphere_pairwise Category Type estimate SE df lower.CL upper.CL t.ratio p.value
## 17 - 200 Left - Right face scrambled 0.1003766 0.2708823 50.83 -0.4434872 0.6442404 0.371 0.7125
##
## Degrees-of-freedom method: satterthwaite
## Confidence level used: 0.95
summary(contra_all_N170_hemi[3], delta = delta_null, level = .9)
summary(contra_all_N170_hemi[3], delta = delta_null, side = ">")
summary(contra_all_N170_hemi[3], delta = delta_null, side = "<")
Comparisons of face-specific N170 amplitude components for scrambled stimuli with different durations:
contra_all_N170_spec[3:4]
## Category_pairwise Duration_pairwise Hemisphere Type estimate SE df lower.CL upper.CL t.ratio p.value
## face - house 17 - 200 Left scrambled 0.15288823 0.2238302 121.85 -0.3550805 0.6608569 0.683 0.9917
## face - house 17 - 200 Right scrambled -0.06344508 0.2238302 121.85 -0.5714138 0.4445236 -0.283 1.0000
##
## Degrees-of-freedom method: satterthwaite
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 2 estimates
## P value adjustment: bonferroni method for 2 tests
summary(contra_all_N170_spec[3:4], delta = delta_null, level = .9)
summary(contra_all_N170_spec[3:4], delta = delta_null, side = ">")
summary(contra_all_N170_spec[3:4], delta = delta_null, side = "<")
Other results
Comparisons of N170 amplitudes between Durations (without multiple comparison corrections):
summary(contra_all_N170[1:8], adjust = "none")
summary(contra_all_N170[1:8], delta = delta_null, level = .9, adjust = "none")
Comparisons of N170 amplitude lateralization between Durations (without multiple comparison corrections):
summary(contra_all_N170_hemi[1:4], adjust = "none")
summary(contra_all_N170_hemi[1:4], delta = delta_null, level = .9)
Comparisons of face-specific N170 amplitudes between Durations (without multiple comparison corrections):
summary(contra_all_N170_spec[1:4], adjust = "none")
summary(contra_all_N170_spec[1:4], delta = delta_null, level = .9, adjust = "none")
Linear mixed model with subjective confidence
Converting response to subjective confidence
df.erp.SC.E2 <- {
df.erp.E2 %>%
filter(Type == "intact") %>% # only keep the intact trials
mutate(Confidence = if_else(Response %in% c("11", "55"), "high",
if_else(Response %in% c("12", "54"), "low",
if_else(substring(Response, 2) == "3", "guess", "NA"))),
DuraConf = paste(Duration, Confidence, sep = "_")) %>%
filter(Confidence != "NA",
!(DuraConf %in% c("200_guess", "200_low"))) %>%
mutate(Confidence = as.factor(Confidence),
DuraConf = as.factor(DuraConf))
}
Set the successive difference coding
# successive difference coding
df.erp.SC.E2 <- sdif_coding_E205_erp(df.erp.SC.E2)
Amplitudes of the P1
# only keep the data for amplitudes of the P1 (correct and incorrect erps)
df.erp.P1.SC.E2 <- {
df.erp.SC.E2 %>%
filter(Component == "P1")
}
if (saveCSV) {
output.erp.P1.SC.E2 <- file.path(folder_nesi, "E205_erp_P1_SC.RData")
save(df.erp.P1.SC.E2, file = output.erp.P1.SC.E2)
}
# the structure of this dataset
head(df.erp.P1.SC.E2, 10)
The maximal model
lmm.max.P1.SC.E2 <- lmer(MeanAmp ~ Category * Hemisphere * DuraConf +
(1 + Cate_C + Hemi_C + Cate_Hemi | SubjCode) +
(1 + Hemi_C | Stimuli),
data = df.erp.P1.SC.E2,
verbose = TRUE,
control = lmerControl(optimizer = "bobyqa")) # nloptwrap Nelder_Mead
## start par. = 1 0 1 1 0 0 0 1 0 0 1 0 1 fn = 221091.6
## At return
## eval: 1576 fn: 220102.61 par: 0.0378072 0.00819983 0.00000 0.338144 -0.0678822 0.0508083 -0.000956374 0.0773424 0.0621299 0.0254038 0.214512 -0.0111326 4.18013e-06
print(summary(lmm.max.P1.SC.E2), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Category * Hemisphere * DuraConf + (1 + Cate_C + Hemi_C + Cate_Hemi | SubjCode) + (1 + Hemi_C | Stimuli)
## Data: df.erp.P1.SC.E2
## Control: lmerControl(optimizer = "bobyqa")
##
## REML criterion at convergence: 220102.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.6020 -0.5787 -0.0155 0.5775 8.4272
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Stimuli (Intercept) 0.021529 0.14673
## Hemi_C 0.001013 0.03182 1.00
## SubjCode (Intercept) 1.722143 1.31230
## Cate_C 0.159498 0.39937 -0.66
## Hemi_C 0.790078 0.88886 0.22 0.06
## Cate_Hemi 0.011600 0.10770 -0.03 0.71 -0.14
## Residual 15.061410 3.88090
## Number of obs: 39600, groups: Stimuli, 80; SubjCode, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.248e+00 2.495e-01 2.730e+01 5.002 2.94e-05 ***
## Category2-1 5.822e-02 9.297e-02 3.683e+01 0.626 0.5350
## Hemisphere2-1 5.166e-02 1.733e-01 2.748e+01 0.298 0.7678
## DuraConf2-1 1.051e-01 6.125e-02 2.983e+04 1.716 0.0862 .
## DuraConf3-2 2.647e-01 6.564e-02 1.776e+04 4.033 5.54e-05 ***
## DuraConf4-3 3.971e-02 6.696e-02 2.288e+04 0.593 0.5532
## Category2-1:Hemisphere2-1 -9.480e-02 8.859e-02 2.180e+02 -1.070 0.2857
## Category2-1:DuraConf2-1 -2.347e-01 1.209e-01 1.212e+04 -1.941 0.0523 .
## Category2-1:DuraConf3-2 3.841e-02 1.263e-01 2.598e+03 0.304 0.7612
## Category2-1:DuraConf4-3 5.656e-02 1.312e-01 6.594e+03 0.431 0.6663
## Hemisphere2-1:DuraConf2-1 1.840e-02 1.205e-01 3.585e+04 0.153 0.8786
## Hemisphere2-1:DuraConf3-2 7.275e-02 1.272e-01 2.436e+04 0.572 0.5673
## Hemisphere2-1:DuraConf4-3 5.993e-02 1.314e-01 2.522e+04 0.456 0.6483
## Category2-1:Hemisphere2-1:DuraConf2-1 1.322e-01 2.334e-01 2.924e+04 0.566 0.5713
## Category2-1:Hemisphere2-1:DuraConf3-2 -1.058e-01 2.330e-01 1.267e+04 -0.454 0.6498
## Category2-1:Hemisphere2-1:DuraConf4-3 -1.807e-01 2.513e-01 2.114e+04 -0.719 0.4721
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## optimizer (bobyqa) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
The zero-correlation-parameter model
lmm.zcp.P1.SC.E2 <- update(lmm.max.P1.SC.E2,
formula = MeanAmp ~ Category * Hemisphere * DuraConf +
(1 + Cate_C + Hemi_C + Cate_Hemi || SubjCode) +
(1 + Hemi_C || Stimuli))
## start par. = 1 1 1 1 1 1 fn = 221091.6
## At return
## eval: 420 fn: 220118.27 par: 2.50536e-07 0.0375916 3.48909e-06 0.229004 0.102779 0.338215
print(summary(lmm.zcp.P1.SC.E2), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Category + Hemisphere + DuraConf + (1 + Cate_C + Hemi_C + Cate_Hemi || SubjCode) + (1 + Hemi_C || Stimuli) + Category:Hemisphere + Category:DuraConf + Hemisphere:DuraConf + Category:Hemisphere:DuraConf
## Data: df.erp.P1.SC.E2
## Control: lmerControl(optimizer = "bobyqa")
##
## REML criterion at convergence: 220118.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.604 -0.578 -0.015 0.578 8.456
##
## Random effects:
## Groups Name Variance Std.Dev.
## Stimuli Hemi_C 0.00000 0.0000
## Stimuli.1 (Intercept) 0.02129 0.1459
## SubjCode Cate_Hemi 0.00000 0.0000
## SubjCode.1 Hemi_C 0.78993 0.8888
## SubjCode.2 Cate_C 0.15911 0.3989
## SubjCode.3 (Intercept) 1.72300 1.3126
## Residual 15.06265 3.8811
## Number of obs: 39600, groups: Stimuli, 80; SubjCode, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.249e+00 2.495e-01 2.730e+01 5.007 2.91e-05 ***
## Category2-1 6.073e-02 9.286e-02 3.677e+01 0.654 0.5172
## Hemisphere2-1 5.095e-02 1.732e-01 2.745e+01 0.294 0.7709
## DuraConf2-1 1.039e-01 6.148e-02 3.568e+04 1.690 0.0911 .
## DuraConf3-2 2.724e-01 6.610e-02 2.887e+04 4.121 3.78e-05 ***
## DuraConf4-3 3.293e-02 6.739e-02 3.365e+04 0.489 0.6251
## Category2-1:Hemisphere2-1 -9.713e-02 8.593e-02 3.917e+04 -1.130 0.2584
## Category2-1:DuraConf2-1 -2.413e-01 1.220e-01 2.111e+04 -1.979 0.0479 *
## Category2-1:DuraConf3-2 6.113e-02 1.290e-01 5.815e+03 0.474 0.6357
## Category2-1:DuraConf4-3 4.158e-02 1.330e-01 1.365e+04 0.313 0.7546
## Hemisphere2-1:DuraConf2-1 2.186e-02 1.205e-01 3.789e+04 0.181 0.8561
## Hemisphere2-1:DuraConf3-2 6.451e-02 1.272e-01 2.990e+04 0.507 0.6121
## Hemisphere2-1:DuraConf4-3 6.541e-02 1.315e-01 3.598e+04 0.497 0.6189
## Category2-1:Hemisphere2-1:DuraConf2-1 1.270e-01 2.332e-01 3.936e+04 0.545 0.5860
## Category2-1:Hemisphere2-1:DuraConf3-2 -1.051e-01 2.325e-01 3.908e+04 -0.452 0.6513
## Category2-1:Hemisphere2-1:DuraConf4-3 -1.768e-01 2.512e-01 3.933e+04 -0.704 0.4814
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## optimizer (bobyqa) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
The reduced model
summary(rePCA(lmm.zcp.P1.SC.E2))
## $Stimuli
## Importance of components:
## [,1] [,2]
## Standard deviation 0.03759 0
## Proportion of Variance 1.00000 0
## Cumulative Proportion 1.00000 1
##
## $SubjCode
## Importance of components:
## [,1] [,2] [,3] [,4]
## Standard deviation 0.3382 0.2290 0.10278 0
## Proportion of Variance 0.6448 0.2956 0.05955 0
## Cumulative Proportion 0.6448 0.9405 1.00000 1
lmm.rdc.P1.SC.E2 <- update(lmm.zcp.P1.SC.E2,
formula = MeanAmp ~ Category * Hemisphere * DuraConf +
(1 + Cate_C + Hemi_C || SubjCode) +
(1 | Stimuli))
## start par. = 1 1 1 1 fn = 220659.7
## At return
## eval: 217 fn: 220118.27 par: 0.0375915 0.229001 0.102778 0.338215
print(summary(lmm.rdc.P1.SC.E2), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Category + Hemisphere + DuraConf + (1 + Cate_C + Hemi_C || SubjCode) + (1 | Stimuli) + Category:Hemisphere + Category:DuraConf + Hemisphere:DuraConf + Category:Hemisphere:DuraConf
## Data: df.erp.P1.SC.E2
## Control: lmerControl(optimizer = "bobyqa")
##
## REML criterion at convergence: 220118.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.604 -0.578 -0.015 0.578 8.456
##
## Random effects:
## Groups Name Variance Std.Dev.
## Stimuli (Intercept) 0.02129 0.1459
## SubjCode Hemi_C 0.78991 0.8888
## SubjCode.1 Cate_C 0.15911 0.3989
## SubjCode.2 (Intercept) 1.72300 1.3126
## Residual 15.06265 3.8811
## Number of obs: 39600, groups: Stimuli, 80; SubjCode, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.249e+00 2.495e-01 2.730e+01 5.007 2.91e-05 ***
## Category2-1 6.073e-02 9.286e-02 3.677e+01 0.654 0.5172
## Hemisphere2-1 5.095e-02 1.732e-01 2.746e+01 0.294 0.7709
## DuraConf2-1 1.039e-01 6.148e-02 3.568e+04 1.690 0.0911 .
## DuraConf3-2 2.724e-01 6.610e-02 2.887e+04 4.121 3.78e-05 ***
## DuraConf4-3 3.293e-02 6.739e-02 3.365e+04 0.489 0.6251
## Category2-1:Hemisphere2-1 -9.713e-02 8.593e-02 3.917e+04 -1.130 0.2584
## Category2-1:DuraConf2-1 -2.413e-01 1.220e-01 2.111e+04 -1.979 0.0479 *
## Category2-1:DuraConf3-2 6.113e-02 1.290e-01 5.815e+03 0.474 0.6357
## Category2-1:DuraConf4-3 4.158e-02 1.330e-01 1.365e+04 0.313 0.7546
## Hemisphere2-1:DuraConf2-1 2.186e-02 1.205e-01 3.789e+04 0.181 0.8561
## Hemisphere2-1:DuraConf3-2 6.451e-02 1.272e-01 2.990e+04 0.507 0.6121
## Hemisphere2-1:DuraConf4-3 6.541e-02 1.315e-01 3.598e+04 0.497 0.6189
## Category2-1:Hemisphere2-1:DuraConf2-1 1.270e-01 2.332e-01 3.936e+04 0.545 0.5860
## Category2-1:Hemisphere2-1:DuraConf3-2 -1.051e-01 2.325e-01 3.908e+04 -0.452 0.6513
## Category2-1:Hemisphere2-1:DuraConf4-3 -1.768e-01 2.512e-01 3.933e+04 -0.704 0.4814
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lmm.rdc.P1.SC.E2, lmm.zcp.P1.SC.E2, refit = FALSE)
The extended model
lmm.etd.P1.SC.E2 <- update(lmm.rdc.P1.SC.E2,
formula = MeanAmp ~ Category * Hemisphere * DuraConf +
(1 + Cate_C + Hemi_C | SubjCode) +
(1 | Stimuli))
## start par. = 1 1 0 0 1 0 1 fn = 220659.7
## At return
## eval: 554 fn: 220104.36 par: 0.0376817 0.338067 -0.0679950 0.0508027 0.0772712 0.0621654 0.214439
print(summary(lmm.etd.P1.SC.E2), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Category + Hemisphere + DuraConf + (1 + Cate_C + Hemi_C | SubjCode) + (1 | Stimuli) + Category:Hemisphere + Category:DuraConf + Hemisphere:DuraConf + Category:Hemisphere:DuraConf
## Data: df.erp.P1.SC.E2
## Control: lmerControl(optimizer = "bobyqa")
##
## REML criterion at convergence: 220104.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.6075 -0.5789 -0.0158 0.5770 8.4361
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Stimuli (Intercept) 0.02139 0.1462
## SubjCode (Intercept) 1.72149 1.3121
## Cate_C 0.15958 0.3995 -0.66
## Hemi_C 0.78972 0.8887 0.22 0.06
## Residual 15.06256 3.8811
## Number of obs: 39600, groups: Stimuli, 80; SubjCode, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.247e+00 2.494e-01 2.730e+01 5.001 2.95e-05 ***
## Category2-1 5.783e-02 9.295e-02 3.674e+01 0.622 0.5377
## Hemisphere2-1 5.142e-02 1.732e-01 2.745e+01 0.297 0.7688
## DuraConf2-1 1.059e-01 6.127e-02 3.012e+04 1.728 0.0840 .
## DuraConf3-2 2.622e-01 6.568e-02 1.817e+04 3.992 6.58e-05 ***
## DuraConf4-3 4.173e-02 6.698e-02 2.308e+04 0.623 0.5333
## Category2-1:Hemisphere2-1 -9.602e-02 8.592e-02 3.911e+04 -1.118 0.2638
## Category2-1:DuraConf2-1 -2.375e-01 1.210e-01 1.240e+04 -1.962 0.0497 *
## Category2-1:DuraConf3-2 3.344e-02 1.265e-01 2.649e+03 0.264 0.7916
## Category2-1:DuraConf4-3 6.136e-02 1.313e-01 6.668e+03 0.467 0.6402
## Hemisphere2-1:DuraConf2-1 2.104e-02 1.205e-01 3.756e+04 0.175 0.8614
## Hemisphere2-1:DuraConf3-2 6.857e-02 1.271e-01 2.857e+04 0.539 0.5896
## Hemisphere2-1:DuraConf4-3 6.261e-02 1.314e-01 3.537e+04 0.476 0.6338
## Category2-1:Hemisphere2-1:DuraConf2-1 1.247e-01 2.332e-01 3.933e+04 0.535 0.5928
## Category2-1:Hemisphere2-1:DuraConf3-2 -9.808e-02 2.324e-01 3.898e+04 -0.422 0.6731
## Category2-1:Hemisphere2-1:DuraConf4-3 -1.824e-01 2.511e-01 3.930e+04 -0.726 0.4677
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lmm.rdc.P1.SC.E2, lmm.etd.P1.SC.E2, refit = FALSE)
The optimal model
lmm.opt.P1.SC.E205 <- lmm.etd.P1.SC.E2
# print(summary(lmm.opt.P1.SC.E2), corr = FALSE)
summary(lmm.opt.P1.SC.E205)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Category + Hemisphere + DuraConf + (1 + Cate_C + Hemi_C | SubjCode) + (1 | Stimuli) + Category:Hemisphere + Category:DuraConf + Hemisphere:DuraConf + Category:Hemisphere:DuraConf
## Data: df.erp.P1.SC.E2
## Control: lmerControl(optimizer = "bobyqa")
##
## REML criterion at convergence: 220104.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.6075 -0.5789 -0.0158 0.5770 8.4361
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Stimuli (Intercept) 0.02139 0.1462
## SubjCode (Intercept) 1.72149 1.3121
## Cate_C 0.15958 0.3995 -0.66
## Hemi_C 0.78972 0.8887 0.22 0.06
## Residual 15.06256 3.8811
## Number of obs: 39600, groups: Stimuli, 80; SubjCode, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.247e+00 2.494e-01 2.730e+01 5.001 2.95e-05 ***
## Category2-1 5.783e-02 9.295e-02 3.674e+01 0.622 0.5377
## Hemisphere2-1 5.142e-02 1.732e-01 2.745e+01 0.297 0.7688
## DuraConf2-1 1.059e-01 6.127e-02 3.012e+04 1.728 0.0840 .
## DuraConf3-2 2.622e-01 6.568e-02 1.817e+04 3.992 6.58e-05 ***
## DuraConf4-3 4.173e-02 6.698e-02 2.308e+04 0.623 0.5333
## Category2-1:Hemisphere2-1 -9.602e-02 8.592e-02 3.911e+04 -1.118 0.2638
## Category2-1:DuraConf2-1 -2.375e-01 1.210e-01 1.240e+04 -1.962 0.0497 *
## Category2-1:DuraConf3-2 3.344e-02 1.265e-01 2.649e+03 0.264 0.7916
## Category2-1:DuraConf4-3 6.136e-02 1.313e-01 6.668e+03 0.467 0.6402
## Hemisphere2-1:DuraConf2-1 2.104e-02 1.205e-01 3.756e+04 0.175 0.8614
## Hemisphere2-1:DuraConf3-2 6.857e-02 1.271e-01 2.857e+04 0.539 0.5896
## Hemisphere2-1:DuraConf4-3 6.261e-02 1.314e-01 3.537e+04 0.476 0.6338
## Category2-1:Hemisphere2-1:DuraConf2-1 1.247e-01 2.332e-01 3.933e+04 0.535 0.5928
## Category2-1:Hemisphere2-1:DuraConf3-2 -9.808e-02 2.324e-01 3.898e+04 -0.422 0.6731
## Category2-1:Hemisphere2-1:DuraConf4-3 -1.824e-01 2.511e-01 3.930e+04 -0.726 0.4677
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
confint(lmm.opt.P1.SC.E205, method = "Wald")
## 2.5 % 97.5 %
## .sig01 NA NA
## .sig02 NA NA
## .sig03 NA NA
## .sig04 NA NA
## .sig05 NA NA
## .sig06 NA NA
## .sig07 NA NA
## .sigma NA NA
## (Intercept) 0.75855980 1.7362675434
## Category2-1 -0.12434173 0.2399972799
## Hemisphere2-1 -0.28803290 0.3908802709
## DuraConf2-1 -0.01421219 0.2259623763
## DuraConf3-2 0.13345670 0.3909005166
## DuraConf4-3 -0.08955719 0.1730168930
## Category2-1:Hemisphere2-1 -0.26442372 0.0723834042
## Category2-1:DuraConf2-1 -0.47465648 -0.0003054915
## Category2-1:DuraConf3-2 -0.21453830 0.2814169651
## Category2-1:DuraConf4-3 -0.19592560 0.3186522135
## Hemisphere2-1:DuraConf2-1 -0.21510192 0.2571773951
## Hemisphere2-1:DuraConf3-2 -0.18060392 0.3177448927
## Hemisphere2-1:DuraConf4-3 -0.19499840 0.3202178645
## Category2-1:Hemisphere2-1:DuraConf2-1 -0.33231528 0.5816902531
## Category2-1:Hemisphere2-1:DuraConf3-2 -0.55367803 0.3575267455
## Category2-1:Hemisphere2-1:DuraConf4-3 -0.67462460 0.3098546508
anova(lmm.opt.P1.SC.E205)
Diagnostic plots
qqplot_lmer(lmm.opt.P1.SC.E205)
Estimated marginal means
emm.P1.SC.E205 <- emmeans(lmm.opt.P1.SC.E205, ~ Category + DuraConf + Hemisphere)
emm.P1.SC.E205 %>%
as_tibble() %>%
arrange(Category, DuraConf)
Plots
# emmip(emm.P1.SC.E205, Category ~ DuraConf | Hemisphere, CIs = TRUE)
duraconf.labels = c("17ms\nguess", "17ms\nlow", "17ms\nhigh", "200ms\nhigh")
p1.conf.LinePlot = {
ggplot(data = data.frame(emm.P1.SC.E205), aes(y = emmean, x = DuraConf, color = Category, group = Category)) + # set the data, varialbes for x and y axises, and the variable for dividing data
geom_point(size = 2) +
geom_line(aes(linetype = Category), size = .85) + # position = "dodge", alpha = .7
scale_linetype_manual(values=c("solid", "dashed")) +
geom_errorbar(mapping = aes(ymin = lower.CL, ymax = upper.CL), linetype = 1, # set the error bar
show.legend = FALSE, width = 0.25, alpha = .5) + # , position = position_dodge(width=0.9)
facet_grid(. ~ Hemisphere, # switch = "x",
labeller = labeller(Hemisphere = hemiLabel)) +
scale_y_reverse(limits = ylimit_p1) + # set the limit for y axis
scale_x_discrete(labels = duraconf.labels) +
# scale_y_continuous(expand= c(0, 0)) + # remove the space between columns and x axis
labs(x = "Duration and Subjective Confidence", y = expression(paste("Amplitudes (", mu, "V)")), fill = "Category") + # set the names for main, x and y axises title = "",
# scale_x_discrete(labels = c("left", "right")) +
# geom_text(label = c("", "", "", "", "***", "", "", ""), color = "red", size = 6, nudge_y = 0.5, nudge_x = 0.5) + # add starts to the significant columns
theme_bw() +
erp_theme
}
# p1.conf.LinePlot.slide <- {
# p1.conf.LinePlot +
# theme(legend.position="right",
# axis.line = element_line(color="black", size = 3))
# }
# ggsave('P1_conf.png', p1.conf.LinePlot.slide, width = 8.5, height = 4)
p1.conf.LinePlot
duration_color <- c(brewer.pal(9, 'Blues')[c(6, 7, 9)], brewer.pal(12, "Paired")[6])
# category_color <- brewer.pal(8, 'Dark2')[c(6,1)]
p1_con_LinePlot = {
data.frame(emm.P1.SC.E205) %>%
mutate(DuraConf = fct_recode(DuraConf,
`33_guess` = "17_guess", `33_low` = "17_low", `33_high` = "17_high", `216_high` = "200_high")) %>%
ggplot(aes(y = emmean, x = DuraConf, group = Category)) + # set the data, varialbes for x and y axises, and the variable for dividing data
geom_point(aes(fill = DuraConf), size=4, shape=21, stroke=0) +
scale_fill_manual(values = duration_color, guide = "none") +
geom_line(aes(linetype = Category)) + # position = "dodge", alpha = .7
scale_linetype_manual(values=c("solid", "dashed")) +
geom_errorbar(mapping = aes(ymin = lower.CL, ymax = upper.CL), linetype = 1, # set the error bar
show.legend = FALSE, width = 0.25, alpha = .5) + # , position = position_dodge(width=0.9)
facet_grid(. ~ Hemisphere, # switch = "x",
labeller = labeller(Hemisphere = hemiLabel)) +
scale_y_reverse(limits = ylimit_p1) + # set the limit for y axis
scale_x_discrete(labels = c("33 ms\nguess", "33 ms\nlow", "33 ms\nhigh", "216 ms\nhigh")) +
# scale_y_continuous(expand= c(0, 0)) + # remove the space between columns and x axis
labs(x = "Duration and Subjective Confidence", y = expression(paste("P1 Amplitudes (", mu, "V)")), fill = "Category") + # set the names for main, x and y axises title = "",
theme_bw() +
erp_theme +
theme(
axis.title.x = element_text(margin = margin(t = 10, r = 0, b = 0, l = 0)),
legend.margin = margin(t = -10),
legend.title=element_text(size=12),
legend.text=element_text(size=12),
)
}
ggsave('P1_con_pub.pdf', p1_con_LinePlot, width = 7, height = 5)
Contrast
Main results
Comparisons of P1 amplitudes for (intact) faces with different durations and confidence (left hemisphere):
contra_SC_P1 <- contrast(emmeans(lmm.opt.P1.SC.E205, ~ DuraConf | Hemisphere + Category), #
"pairwise", adjust = "Bonferroni", infer = TRUE)
contra_SC_P1[1:6]
## contrast Hemisphere Category estimate SE df lower.CL upper.CL t.ratio p.value
## 17_guess - 17_low Left face -0.2452686 0.1354348 37115.10 -0.6025995 0.1120623 -1.811 0.4209
## 17_guess - 17_high Left face -0.4319234 0.1376417 20687.49 -0.7950925 -0.0687543 -3.138 0.0102
## 17_guess - 200_high Left face -0.3660705 0.1425063 36491.93 -0.7420593 0.0099183 -2.569 0.0613
## 17_low - 17_high Left face -0.1866548 0.1094711 17691.79 -0.4755002 0.1021907 -1.705 0.5292
## 17_low - 200_high Left face -0.1208019 0.1171315 37690.70 -0.4298413 0.1882376 -1.031 1.0000
## 17_high - 200_high Left face 0.0658529 0.1097015 36231.43 -0.2235837 0.3552895 0.600 1.0000
##
## Degrees-of-freedom method: satterthwaite
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 6 estimates
## P value adjustment: bonferroni method for 6 tests
summary(contra_SC_P1[1:6], delta = delta_null, level = .9)
summary(contra_SC_P1[1:6], delta = delta_null, side = ">")
summary(contra_SC_P1[1:6], delta = delta_null, side = "<")
Comparisons of P1 amplitudes for (intact) faces with different durations and confidence (right hemisphere):
contra_SC_P1[7:12]
## contrast Hemisphere Category estimate SE df lower.CL upper.CL t.ratio p.value
## 17_guess - 17_low Right face -0.2039626 0.1352749 35346.56 -0.5608726 0.15294740 -1.508 0.7897
## 17_guess - 17_high Right face -0.5082257 0.1370905 15779.36 -0.8699512 -0.14650017 -3.707 0.0013
## 17_guess - 200_high Right face -0.5961750 0.1423218 34614.65 -0.9716781 -0.22067186 -4.189 0.0002
## 17_low - 17_high Right face -0.3042631 0.1089672 13017.88 -0.5917906 -0.01673559 -2.792 0.0315
## 17_low - 200_high Right face -0.3922124 0.1170154 36493.10 -0.7009459 -0.08347888 -3.352 0.0048
## 17_high - 200_high Right face -0.0879493 0.1095568 34348.22 -0.3770051 0.20110645 -0.803 1.0000
##
## Degrees-of-freedom method: satterthwaite
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 6 estimates
## P value adjustment: bonferroni method for 6 tests
summary(contra_SC_P1[7:12], delta = delta_null, level = .9)
summary(contra_SC_P1[7:12], delta = delta_null, side = ">")
summary(contra_SC_P1[7:12], delta = delta_null, side = "<")
Comparisons of face-specific P1 amplitude components for different durations and confidence (left hemisphere):
contra_SC_P1_spec <- contrast(emmeans(lmm.opt.P1.SC.E205, ~ Category + DuraConf | Hemisphere), #
interaction = "pairwise", adjust = "Bonferroni", infer = TRUE)
contra_SC_P1_spec[1:6]
## Category_pairwise DuraConf_pairwise Hemisphere estimate SE df lower.CL upper.CL t.ratio p.value
## face - house 17_guess - 17_low Left -0.29982473 0.1679815 24749.16 -0.7430388 0.1433893 -1.785 0.4458
## face - house 17_guess - 17_high Left -0.21734758 0.1920036 5076.76 -0.7241011 0.2894060 -1.132 1.0000
## face - house 17_guess - 200_high Left -0.06479179 0.1800928 31977.11 -0.5399525 0.4103690 -0.360 1.0000
## face - house 17_low - 17_high Left 0.08247715 0.1716759 7492.37 -0.3705685 0.5355228 0.480 1.0000
## face - house 17_low - 200_high Left 0.23503294 0.1617098 38100.49 -0.1916213 0.6616872 1.453 0.8767
## face - house 17_high - 200_high Left 0.15255579 0.1814463 15301.12 -0.3262084 0.6313200 0.841 1.0000
##
## Degrees-of-freedom method: satterthwaite
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 6 estimates
## P value adjustment: bonferroni method for 6 tests
summary(contra_SC_P1_spec[1:6], delta = delta_null, level = .9)
summary(contra_SC_P1_spec[1:6], delta = delta_null, side = ">")
summary(contra_SC_P1_spec[1:6], delta = delta_null, side = "<")
Comparisons of face-specific P1 amplitude components for different durations and confidence (right hemisphere):
contra_SC_P1_spec[7:12]
## Category_pairwise DuraConf_pairwise Hemisphere estimate SE df lower.CL upper.CL t.ratio p.value
## face - house 17_guess - 17_low Right -0.17513724 0.1680863 25502.49 -0.6186267 0.2683523 -1.042 1.0000
## face - house 17_guess - 17_high Right -0.19073573 0.1923894 5541.46 -0.6984908 0.3170193 -0.991 1.0000
## face - house 17_guess - 200_high Right -0.22056491 0.1800024 31001.27 -0.6954880 0.2543582 -1.225 1.0000
## face - house 17_low - 17_high Right -0.01559849 0.1719311 8082.09 -0.4693087 0.4381117 -0.091 1.0000
## face - house 17_low - 200_high Right -0.04542767 0.1616989 37871.29 -0.4720533 0.3811979 -0.281 1.0000
## face - house 17_high - 200_high Right -0.02982918 0.1818778 17940.88 -0.5097230 0.4500646 -0.164 1.0000
##
## Degrees-of-freedom method: satterthwaite
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 6 estimates
## P value adjustment: bonferroni method for 6 tests
summary(contra_SC_P1_spec[7:12], delta = delta_null, level = .9)
summary(contra_SC_P1_spec[7:12], delta = delta_null, side = ">")
summary(contra_SC_P1_spec[7:12], delta = delta_null, side = "<")
Other results
Comparisons of P1 amplitudes for (intact) faces with different durations and confidence (without multiple comparison corrections):
summary(contra_SC_P1[1:24], adjust = "none")
summary(contra_SC_P1[1:24], adjust = "none", delta = delta_null, level = .9)
Ampitudes of the N170
# only keep the data for amplitudes of the N170 (all erps)
df.erp.N170.SC.E2 <- {
df.erp.SC.E2 %>%
filter(Component == "N170")
}
if (saveCSV) {
output.erp.N170.SC.E2 <- file.path(folder_nesi, "E205_erp_N170_SC.RData")
save(df.erp.N170.SC.E2, file = output.erp.N170.SC.E2)
}
# the structure of this dataset
head(df.erp.N170.SC.E2, 10)
The maximal model
lmm.max.N170.SC.E2 <- lmer(MeanAmp ~ Category * Hemisphere * DuraConf +
(1 + Cate_C + Hemi_C + Cate_Hemi | SubjCode) +
(1 + Hemi_C | Stimuli),
data = df.erp.N170.SC.E2,
verbose = TRUE,
control = lmerControl(optimizer = "bobyqa")) # nloptwrap Nelder_Mead
## start par. = 1 0 1 1 0 0 0 1 0 0 1 0 1 fn = 222383.5
## At return
## eval: 1455 fn: 221500.89 par: 0.0377452 0.00592604 1.27875e-07 0.468916 -0.124988 -0.00676600 -0.0181415 0.172582 -0.0131822 -0.0558530 0.351847 -0.213334 0.184849
print(summary(lmm.max.N170.SC.E2), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Category * Hemisphere * DuraConf + (1 + Cate_C + Hemi_C + Cate_Hemi | SubjCode) + (1 + Hemi_C | Stimuli)
## Data: df.erp.N170.SC.E2
## Control: lmerControl(optimizer = "bobyqa")
##
## REML criterion at convergence: 221500.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.3165 -0.5946 0.0000 0.5895 6.7424
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Stimuli (Intercept) 2.217e-02 0.14888
## Hemi_C 5.464e-04 0.02337 1.00
## SubjCode (Intercept) 3.421e+00 1.84956
## Cate_C 7.064e-01 0.84049 -0.59
## Hemi_C 1.929e+00 1.38903 -0.02 -0.02
## Cate_Hemi 1.293e+00 1.13724 -0.06 -0.12 -0.73
## Residual 1.556e+01 3.94434
## Number of obs: 39600, groups: Stimuli, 80; SubjCode, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.996e+00 3.506e-01 2.714e+01 -5.693 4.68e-06 ***
## Category2-1 1.217e+00 1.682e-01 2.946e+01 7.236 5.24e-08 ***
## Hemisphere2-1 -1.775e-01 2.661e-01 2.726e+01 -0.667 0.51039
## DuraConf2-1 -4.935e-01 6.263e-02 3.817e+04 -7.880 3.35e-15 ***
## DuraConf3-2 -1.805e-01 6.747e-02 3.568e+04 -2.675 0.00747 **
## DuraConf4-3 3.715e-01 6.865e-02 3.638e+04 5.412 6.29e-08 ***
## Category2-1:Hemisphere2-1 3.189e-01 2.322e-01 2.856e+01 1.373 0.18036
## Category2-1:DuraConf2-1 5.042e-01 1.248e-01 3.307e+04 4.039 5.37e-05 ***
## Category2-1:DuraConf3-2 3.517e-01 1.335e-01 1.716e+04 2.635 0.00843 **
## Category2-1:DuraConf4-3 1.770e+00 1.365e-01 2.594e+04 12.971 < 2e-16 ***
## Hemisphere2-1:DuraConf2-1 -1.560e-03 1.242e-01 3.044e+04 -0.013 0.98998
## Hemisphere2-1:DuraConf3-2 -2.377e-01 1.329e-01 1.699e+04 -1.789 0.07358 .
## Hemisphere2-1:DuraConf4-3 -1.685e-01 1.354e-01 1.950e+04 -1.244 0.21353
## Category2-1:Hemisphere2-1:DuraConf2-1 2.736e-01 2.465e-01 1.737e+04 1.110 0.26707
## Category2-1:Hemisphere2-1:DuraConf3-2 -5.786e-02 2.597e-01 4.298e+03 -0.223 0.82372
## Category2-1:Hemisphere2-1:DuraConf4-3 2.674e-01 2.675e-01 8.591e+03 1.000 0.31754
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## optimizer (bobyqa) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
The zero-correlation-parameter model
lmm.zcp.N170.SC.E2 <- update(lmm.max.N170.SC.E2,
formula = MeanAmp ~ Category * Hemisphere * DuraConf +
(1 + Cate_C + Hemi_C + Cate_Hemi || SubjCode) +
(1 + Hemi_C || Stimuli))
## start par. = 1 1 1 1 1 1 fn = 222383.5
## At return
## eval: 464 fn: 221530.50 par: 0.00000 0.0376425 0.288392 0.352220 0.212488 0.469156
print(summary(lmm.zcp.N170.SC.E2), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Category + Hemisphere + DuraConf + (1 + Cate_C + Hemi_C + Cate_Hemi || SubjCode) + (1 + Hemi_C || Stimuli) + Category:Hemisphere + Category:DuraConf + Hemisphere:DuraConf + Category:Hemisphere:DuraConf
## Data: df.erp.N170.SC.E2
## Control: lmerControl(optimizer = "bobyqa")
##
## REML criterion at convergence: 221530.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.3076 -0.5935 -0.0001 0.5902 6.7081
##
## Random effects:
## Groups Name Variance Std.Dev.
## Stimuli Hemi_C 0.00000 0.0000
## Stimuli.1 (Intercept) 0.02205 0.1485
## SubjCode Cate_Hemi 1.29397 1.1375
## SubjCode.1 Hemi_C 1.93012 1.3893
## SubjCode.2 Cate_C 0.70247 0.8381
## SubjCode.3 (Intercept) 3.42446 1.8505
## Residual 15.55810 3.9444
## Number of obs: 39600, groups: Stimuli, 80; SubjCode, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.995e+00 3.508e-01 2.714e+01 -5.687 4.76e-06 ***
## Category2-1 1.217e+00 1.677e-01 2.946e+01 7.253 5.00e-08 ***
## Hemisphere2-1 -1.757e-01 2.662e-01 2.726e+01 -0.660 0.51480
## DuraConf2-1 -4.958e-01 6.270e-02 3.886e+04 -7.909 2.67e-15 ***
## DuraConf3-2 -1.785e-01 6.761e-02 3.821e+04 -2.640 0.00829 **
## DuraConf4-3 3.698e-01 6.879e-02 3.855e+04 5.376 7.66e-08 ***
## Category2-1:Hemisphere2-1 3.189e-01 2.323e-01 2.856e+01 1.373 0.18044
## Category2-1:DuraConf2-1 5.018e-01 1.251e-01 3.672e+04 4.011 6.05e-05 ***
## Category2-1:DuraConf3-2 3.702e-01 1.343e-01 2.615e+04 2.757 0.00584 **
## Category2-1:DuraConf4-3 1.759e+00 1.371e-01 3.347e+04 12.835 < 2e-16 ***
## Hemisphere2-1:DuraConf2-1 3.713e-03 1.248e-01 3.817e+04 0.030 0.97626
## Hemisphere2-1:DuraConf3-2 -2.309e-01 1.342e-01 3.377e+04 -1.721 0.08521 .
## Hemisphere2-1:DuraConf4-3 -1.775e-01 1.368e-01 3.704e+04 -1.298 0.19430
## Category2-1:Hemisphere2-1:DuraConf2-1 2.960e-01 2.486e-01 3.101e+04 1.191 0.23377
## Category2-1:Hemisphere2-1:DuraConf3-2 -3.908e-02 2.652e-01 1.341e+04 -0.147 0.88287
## Category2-1:Hemisphere2-1:DuraConf4-3 2.456e-01 2.718e-01 2.410e+04 0.904 0.36614
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## optimizer (bobyqa) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
anova(lmm.max.N170.SC.E2, lmm.zcp.N170.SC.E2, refit = FALSE)
The reduced model
summary(rePCA(lmm.zcp.N170.SC.E2))
## $Stimuli
## Importance of components:
## [,1] [,2]
## Standard deviation 0.03764 0
## Proportion of Variance 1.00000 0
## Cumulative Proportion 1.00000 1
##
## $SubjCode
## Importance of components:
## [,1] [,2] [,3] [,4]
## Standard deviation 0.4692 0.3522 0.2884 0.21249
## Proportion of Variance 0.4658 0.2626 0.1760 0.09556
## Cumulative Proportion 0.4658 0.7284 0.9044 1.00000
lmm.rdc.N170.SC.E2 <- update(lmm.zcp.N170.SC.E2,
formula = MeanAmp ~ Category * Hemisphere * DuraConf +
(1 + Cate_C + Hemi_C + Cate_Hemi || SubjCode) +
(1 | Stimuli))
## start par. = 1 1 1 1 1 fn = 222042.5
## At return
## eval: 394 fn: 221530.50 par: 0.0376426 0.288391 0.352220 0.212489 0.469156
print(summary(lmm.rdc.N170.SC.E2), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Category + Hemisphere + DuraConf + (1 + Cate_C + Hemi_C + Cate_Hemi || SubjCode) + (1 | Stimuli) + Category:Hemisphere + Category:DuraConf + Hemisphere:DuraConf + Category:Hemisphere:DuraConf
## Data: df.erp.N170.SC.E2
## Control: lmerControl(optimizer = "bobyqa")
##
## REML criterion at convergence: 221530.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.3076 -0.5935 -0.0001 0.5902 6.7081
##
## Random effects:
## Groups Name Variance Std.Dev.
## Stimuli (Intercept) 0.02205 0.1485
## SubjCode Cate_Hemi 1.29396 1.1375
## SubjCode.1 Hemi_C 1.93012 1.3893
## SubjCode.2 Cate_C 0.70247 0.8381
## SubjCode.3 (Intercept) 3.42446 1.8505
## Residual 15.55810 3.9444
## Number of obs: 39600, groups: Stimuli, 80; SubjCode, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.995e+00 3.508e-01 2.714e+01 -5.687 4.76e-06 ***
## Category2-1 1.217e+00 1.677e-01 2.946e+01 7.253 5.00e-08 ***
## Hemisphere2-1 -1.757e-01 2.662e-01 2.726e+01 -0.660 0.51480
## DuraConf2-1 -4.958e-01 6.270e-02 3.886e+04 -7.909 2.67e-15 ***
## DuraConf3-2 -1.785e-01 6.761e-02 3.821e+04 -2.640 0.00829 **
## DuraConf4-3 3.698e-01 6.879e-02 3.855e+04 5.376 7.66e-08 ***
## Category2-1:Hemisphere2-1 3.189e-01 2.323e-01 2.856e+01 1.373 0.18044
## Category2-1:DuraConf2-1 5.018e-01 1.251e-01 3.672e+04 4.011 6.05e-05 ***
## Category2-1:DuraConf3-2 3.702e-01 1.343e-01 2.615e+04 2.757 0.00584 **
## Category2-1:DuraConf4-3 1.759e+00 1.371e-01 3.347e+04 12.835 < 2e-16 ***
## Hemisphere2-1:DuraConf2-1 3.713e-03 1.248e-01 3.817e+04 0.030 0.97626
## Hemisphere2-1:DuraConf3-2 -2.309e-01 1.342e-01 3.377e+04 -1.721 0.08521 .
## Hemisphere2-1:DuraConf4-3 -1.775e-01 1.368e-01 3.704e+04 -1.298 0.19430
## Category2-1:Hemisphere2-1:DuraConf2-1 2.960e-01 2.486e-01 3.101e+04 1.191 0.23377
## Category2-1:Hemisphere2-1:DuraConf3-2 -3.908e-02 2.652e-01 1.341e+04 -0.147 0.88287
## Category2-1:Hemisphere2-1:DuraConf4-3 2.456e-01 2.718e-01 2.410e+04 0.904 0.36614
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lmm.rdc.N170.SC.E2, lmm.zcp.N170.SC.E2, refit = FALSE)
The extended model
lmm.etd.N170.SC.E2 <- update(lmm.rdc.N170.SC.E2,
formula = MeanAmp ~ Category * Hemisphere * DuraConf +
(1 + Cate_C + Hemi_C + Cate_Hemi | SubjCode) +
(1 | Stimuli))
## start par. = 1 1 0 0 0 1 0 0 1 0 1 fn = 222042.5
## At return
## eval: 1146 fn: 221501.14 par: 0.0376527 0.468917 -0.124977 -0.00678833 -0.0181633 0.172587 -0.0131156 -0.0558089 0.351828 -0.213357 0.184805
print(summary(lmm.rdc.N170.SC.E2), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Category + Hemisphere + DuraConf + (1 + Cate_C + Hemi_C + Cate_Hemi || SubjCode) + (1 | Stimuli) + Category:Hemisphere + Category:DuraConf + Hemisphere:DuraConf + Category:Hemisphere:DuraConf
## Data: df.erp.N170.SC.E2
## Control: lmerControl(optimizer = "bobyqa")
##
## REML criterion at convergence: 221530.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.3076 -0.5935 -0.0001 0.5902 6.7081
##
## Random effects:
## Groups Name Variance Std.Dev.
## Stimuli (Intercept) 0.02205 0.1485
## SubjCode Cate_Hemi 1.29396 1.1375
## SubjCode.1 Hemi_C 1.93012 1.3893
## SubjCode.2 Cate_C 0.70247 0.8381
## SubjCode.3 (Intercept) 3.42446 1.8505
## Residual 15.55810 3.9444
## Number of obs: 39600, groups: Stimuli, 80; SubjCode, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.995e+00 3.508e-01 2.714e+01 -5.687 4.76e-06 ***
## Category2-1 1.217e+00 1.677e-01 2.946e+01 7.253 5.00e-08 ***
## Hemisphere2-1 -1.757e-01 2.662e-01 2.726e+01 -0.660 0.51480
## DuraConf2-1 -4.958e-01 6.270e-02 3.886e+04 -7.909 2.67e-15 ***
## DuraConf3-2 -1.785e-01 6.761e-02 3.821e+04 -2.640 0.00829 **
## DuraConf4-3 3.698e-01 6.879e-02 3.855e+04 5.376 7.66e-08 ***
## Category2-1:Hemisphere2-1 3.189e-01 2.323e-01 2.856e+01 1.373 0.18044
## Category2-1:DuraConf2-1 5.018e-01 1.251e-01 3.672e+04 4.011 6.05e-05 ***
## Category2-1:DuraConf3-2 3.702e-01 1.343e-01 2.615e+04 2.757 0.00584 **
## Category2-1:DuraConf4-3 1.759e+00 1.371e-01 3.347e+04 12.835 < 2e-16 ***
## Hemisphere2-1:DuraConf2-1 3.713e-03 1.248e-01 3.817e+04 0.030 0.97626
## Hemisphere2-1:DuraConf3-2 -2.309e-01 1.342e-01 3.377e+04 -1.721 0.08521 .
## Hemisphere2-1:DuraConf4-3 -1.775e-01 1.368e-01 3.704e+04 -1.298 0.19430
## Category2-1:Hemisphere2-1:DuraConf2-1 2.960e-01 2.486e-01 3.101e+04 1.191 0.23377
## Category2-1:Hemisphere2-1:DuraConf3-2 -3.908e-02 2.652e-01 1.341e+04 -0.147 0.88287
## Category2-1:Hemisphere2-1:DuraConf4-3 2.456e-01 2.718e-01 2.410e+04 0.904 0.36614
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lmm.rdc.N170.SC.E2, lmm.etd.N170.SC.E2, refit = FALSE)
The optimal model
lmm.opt.N170.SC.E205 <- lmm.etd.N170.SC.E2
print(summary(lmm.opt.N170.SC.E205), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Category + Hemisphere + DuraConf + (1 + Cate_C + Hemi_C + Cate_Hemi | SubjCode) + (1 | Stimuli) + Category:Hemisphere + Category:DuraConf + Hemisphere:DuraConf + Category:Hemisphere:DuraConf
## Data: df.erp.N170.SC.E2
## Control: lmerControl(optimizer = "bobyqa")
##
## REML criterion at convergence: 221501.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.3165 -0.5946 -0.0004 0.5897 6.7418
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Stimuli (Intercept) 0.02206 0.1485
## SubjCode (Intercept) 3.42096 1.8496
## Cate_C 0.70642 0.8405 -0.59
## Hemi_C 1.92921 1.3890 -0.02 -0.02
## Cate_Hemi 1.29317 1.1372 -0.06 -0.12 -0.73
## Residual 15.55805 3.9444
## Number of obs: 39600, groups: Stimuli, 80; SubjCode, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.996e+00 3.506e-01 2.714e+01 -5.693 4.68e-06 ***
## Category2-1 1.217e+00 1.681e-01 2.945e+01 7.236 5.24e-08 ***
## Hemisphere2-1 -1.775e-01 2.661e-01 2.725e+01 -0.667 0.51030
## DuraConf2-1 -4.936e-01 6.263e-02 3.817e+04 -7.882 3.31e-15 ***
## DuraConf3-2 -1.805e-01 6.747e-02 3.568e+04 -2.674 0.00749 **
## DuraConf4-3 3.715e-01 6.865e-02 3.638e+04 5.411 6.29e-08 ***
## Category2-1:Hemisphere2-1 3.184e-01 2.321e-01 2.853e+01 1.372 0.18085
## Category2-1:DuraConf2-1 5.042e-01 1.248e-01 3.307e+04 4.040 5.36e-05 ***
## Category2-1:DuraConf3-2 3.517e-01 1.335e-01 1.717e+04 2.634 0.00844 **
## Category2-1:DuraConf4-3 1.770e+00 1.365e-01 2.596e+04 12.970 < 2e-16 ***
## Hemisphere2-1:DuraConf2-1 -3.097e-03 1.242e-01 3.044e+04 -0.025 0.98010
## Hemisphere2-1:DuraConf3-2 -2.384e-01 1.329e-01 1.699e+04 -1.794 0.07281 .
## Hemisphere2-1:DuraConf4-3 -1.673e-01 1.354e-01 1.950e+04 -1.235 0.21683
## Category2-1:Hemisphere2-1:DuraConf2-1 2.733e-01 2.465e-01 1.737e+04 1.109 0.26762
## Category2-1:Hemisphere2-1:DuraConf3-2 -5.872e-02 2.597e-01 4.297e+03 -0.226 0.82116
## Category2-1:Hemisphere2-1:DuraConf4-3 2.686e-01 2.675e-01 8.591e+03 1.004 0.31538
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lmm.opt.N170.SC.E205)
Diagnostic plots
qqplot_lmer(lmm.opt.N170.SC.E205)
Estimated marginal means
emm.N170.SC.E205 <- emmeans(lmm.opt.N170.SC.E205, ~ Category + DuraConf + Hemisphere)
emm.N170.SC.E205 %>%
as_tibble() %>%
arrange(Category, DuraConf)
Plots
# emmip(emm.N170.SC.E205, Category ~ DuraConf | Hemisphere, CIs = TRUE) +
# scale_y_continuous(trans = "reverse")
duraconf.labels = c("17ms\nguess", "17ms\nlow", "17ms\nhigh", "200ms\nhigh")
n170.conf.LinePlot = {
ggplot(data = data.frame(emm.N170.SC.E205), aes(y = emmean, x = DuraConf, color = Category, group = Category)) + # set the data, varialbes for x and y axises, and the variable for dividing data
geom_point(size = 2) +
geom_line(aes(linetype = Category), size = .85) + # position = "dodge", alpha = .7
scale_linetype_manual(values=c("solid", "dashed")) +
geom_errorbar(mapping = aes(ymin = lower.CL, ymax = upper.CL), linetype = 1, # set the error bar
show.legend = FALSE, width = 0.25, alpha = .5) + # , position = position_dodge(width=0.9)
facet_grid(. ~ Hemisphere, # switch = "x",
labeller = labeller(Hemisphere = hemiLabel)) +
scale_y_reverse(limits = ylimit_n170, breaks = seq(1, -5, -1)) + # set the limit for y axis
scale_x_discrete(labels = duraconf.labels) +
# scale_y_continuous(expand= c(0, 0)) + # remove the space between columns and x axis
labs(x = "Duration and Subjective Confidence", y = expression(paste("Amplitudes (", mu, "V)")), fill = "Category") + # set the names for main, x and y axises title = "",
# scale_x_discrete(labels = c("left", "right")) +
# geom_text(label = c("", "", "", "", "***", "", "", ""), color = "red", size = 6, nudge_y = 0.5, nudge_x = 0.5) + # add starts to the significant columns
theme_bw() +
erp_theme
}
# n170.conf.LinePlot.slide <- {
# n170.conf.LinePlot +
# scale_y_reverse(limits = c(0.8, -4.8), breaks = seq(0, -4, -1)) + # set the limit for y axis
# theme(legend.position="right",
# axis.line = element_line(color="black", size = 3))
# }
# ggsave('N170_confi.png', n170.conf.LinePlot.slide, width = 8.5, height = 4)
n170.conf.LinePlot
duration_color <- c(brewer.pal(9, 'Blues')[c(6, 7, 9)], brewer.pal(12, "Paired")[6])
# category_color <- brewer.pal(8, 'Dark2')[c(6,1)]
n170_con_LinePlot = {
data.frame(emm.N170.SC.E205) %>%
mutate(DuraConf = fct_recode(DuraConf,
`33_guess` = "17_guess", `33_low` = "17_low", `33_high` = "17_high", `216_high` = "200_high")) %>%
ggplot(aes(y = emmean, x = DuraConf, group = Category)) + # set the data, varialbes for x and y axises, and the variable for dividing data
geom_point(aes(fill = DuraConf), size=4, shape=21, stroke=0) +
scale_fill_manual(values = duration_color, guide = "none") +
geom_line(aes(linetype = Category)) + # position = "dodge", alpha = .7
scale_linetype_manual(values=c("solid", "dashed")) +
geom_errorbar(mapping = aes(ymin = lower.CL, ymax = upper.CL), linetype = 1, # set the error bar
show.legend = FALSE, width = 0.25, alpha = .5) + # , position = position_dodge(width=0.9)
facet_grid(. ~ Hemisphere, # switch = "x",
labeller = labeller(Hemisphere = hemiLabel)) +
scale_y_reverse(limits = ylimit_n170, breaks = seq(1, -5, -1)) + # set the limit for y axis
scale_x_discrete(labels = c("33 ms\nguess", "33 ms\nlow", "33 ms\nhigh", "216 ms\nhigh")) +
# scale_y_continuous(expand= c(0, 0)) + # remove the space between columns and x axis
labs(x = "Duration and Subjective Confidence", y = expression(paste("N170 Amplitudes (", mu, "V)")), fill = "Category") + # set the names for main, x and y axises title = "",
theme_bw() +
erp_theme +
theme(
axis.title.x = element_text(margin = margin(t = 10, r = 0, b = 0, l = 0)),
legend.margin = margin(t = -10),
legend.title=element_text(size=12),
legend.text=element_text(size=12),
)
}
ggsave('N170_con_pub.pdf', n170_con_LinePlot, width = 7, height = 5)
Contrast
Main results
Comparisons of N170 amplitudes for (intact) faces with different durations and confidence (left hemisphere):
contra_SC_N170 <- contrast(emmeans(lmm.opt.N170.SC.E205, ~ DuraConf | Hemisphere + Category), #
"pairwise", adjust = "Bonferroni", infer = TRUE)
contra_SC_N170[1:6]
## contrast Hemisphere Category estimate SE df lower.CL upper.CL t.ratio p.value
## 17_guess - 17_low Left face 0.6758619 0.1395425 38713.15 0.3076941 1.0440298 4.843 <.0001
## 17_guess - 17_high Left face 0.9276422 0.1460076 29514.59 0.5424106 1.3128739 6.353 <.0001
## 17_guess - 200_high Left face 1.2904477 0.1469225 38302.54 0.9028082 1.6780872 8.783 <.0001
## 17_low - 17_high Left face 0.2517803 0.1169049 27206.64 -0.0566675 0.5602280 2.154 0.1876
## 17_low - 200_high Left face 0.6145858 0.1203932 38772.85 0.2969414 0.9322302 5.105 <.0001
## 17_high - 200_high Left face 0.3628055 0.1131308 38093.91 0.0643219 0.6612891 3.207 0.0081
##
## Degrees-of-freedom method: satterthwaite
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 6 estimates
## P value adjustment: bonferroni method for 6 tests
summary(contra_SC_N170[1:6], delta = delta_null, level = .9)
summary(contra_SC_N170[1:6], delta = delta_null, side = ">")
summary(contra_SC_N170[1:6], delta = delta_null, side = "<")
Comparisons of N170 amplitudes for (intact) faces with different durations and confidence (right hemisphere):
contra_SC_N170[7:12]
## contrast Hemisphere Category estimate SE df lower.CL upper.CL t.ratio p.value
## 17_guess - 17_low Right face 0.8155917 0.1393891 37587.07 0.4478279 1.1833556 5.851 <.0001
## 17_guess - 17_high Right face 1.2764073 0.1454449 22943.21 0.8926529 1.6601616 8.776 <.0001
## 17_guess - 200_high Right face 1.9407591 0.1467338 36887.28 1.5536167 2.3279014 13.226 <.0001
## 17_low - 17_high Right face 0.4608155 0.1163863 20095.74 0.1537282 0.7679028 3.959 0.0005
## 17_low - 200_high Right face 1.1251673 0.1202709 37878.28 0.8078451 1.4424896 9.355 <.0001
## 17_high - 200_high Right face 0.6643518 0.1129761 36540.10 0.3662755 0.9624281 5.880 <.0001
##
## Degrees-of-freedom method: satterthwaite
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 6 estimates
## P value adjustment: bonferroni method for 6 tests
summary(contra_SC_N170[7:12], delta = delta_null, level = .9)
summary(contra_SC_N170[7:12], delta = delta_null, side = ">")
summary(contra_SC_N170[7:12], delta = delta_null, side = "<")
Comparisons of N170 amplitude lateralization for (intact) faces with different durations and confidence:
contra_SC_N170_hemi <- contrast(emmeans(lmm.opt.N170.SC.E205, ~ DuraConf + Hemisphere | Category), #
interaction = "pairwise", adjust = "Bonferroni", infer = TRUE)
contra_SC_N170_hemi[1:6]
## DuraConf_pairwise Hemisphere_pairwise Category estimate SE df lower.CL upper.CL t.ratio p.value
## 17_guess - 17_low Left - Right face -0.1397298 0.1969420 37121.29 -0.6593412 0.3798816 -0.709 1.0000
## 17_guess - 17_high Left - Right face -0.3487650 0.2052731 21129.07 -0.8903792 0.1928491 -1.699 0.5360
## 17_guess - 200_high Left - Right face -0.6503113 0.2073552 36333.86 -1.1973977 -0.1032250 -3.136 0.0103
## 17_low - 17_high Left - Right face -0.2090352 0.1642874 18263.34 -0.6425150 0.2244445 -1.272 1.0000
## 17_low - 200_high Left - Right face -0.5105815 0.1700071 37522.46 -0.9591277 -0.0620354 -3.003 0.0160
## 17_high - 200_high Left - Right face -0.3015463 0.1596700 35954.30 -0.7228200 0.1197275 -1.889 0.3538
##
## Degrees-of-freedom method: satterthwaite
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 6 estimates
## P value adjustment: bonferroni method for 6 tests
summary(contra_SC_N170_hemi[1:6], delta = delta_null, level = .9)
summary(contra_SC_N170_hemi[1:6], delta = delta_null, side = ">")
summary(contra_SC_N170_hemi[1:6], delta = delta_null, side = "<")
Comparisons of face-specific N170 amplitude components for different durations and confidence (left hemisphere):
contra_SC_N170_spec <- contrast(emmeans(lmm.opt.N170.SC.E205, ~ Category + DuraConf | Hemisphere), #
interaction = "pairwise", adjust = "Bonferroni", infer = TRUE)
contra_SC_N170_spec[1:6]
## Category_pairwise DuraConf_pairwise Hemisphere estimate SE df lower.CL upper.CL t.ratio p.value
## face - house 17_guess - 17_low Left 0.3676071 0.1758644 29292.96 -0.0963998 0.8316141 2.090 0.2196
## face - house 17_guess - 17_high Left 0.7486509 0.2128594 8673.02 0.1869443 1.3103575 3.517 0.0026
## face - house 17_guess - 200_high Left 2.3845078 0.1868966 34616.71 1.8913980 2.8776175 12.758 <.0001
## face - house 17_low - 17_high Left 0.3810437 0.1874597 11931.55 -0.1136056 0.8756931 2.033 0.2526
## face - house 17_low - 200_high Left 2.0169006 0.1659301 38713.14 1.5791118 2.4546895 12.155 <.0001
## face - house 17_high - 200_high Left 1.6358569 0.1919169 19722.82 1.1294796 2.1422342 8.524 <.0001
##
## Degrees-of-freedom method: satterthwaite
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 6 estimates
## P value adjustment: bonferroni method for 6 tests
summary(contra_SC_N170_spec[1:6], delta = delta_null, level = .9)
summary(contra_SC_N170_spec[1:6], delta = delta_null, side = ">")
summary(contra_SC_N170_spec[1:6], delta = delta_null, side = "<")
Comparisons of face-specific N170 amplitude components for different durations and confidence (right hemisphere):
contra_SC_N170_spec[7:12]
## Category_pairwise DuraConf_pairwise Hemisphere estimate SE df lower.CL upper.CL t.ratio p.value
## face - house 17_guess - 17_low Right 0.6408724 0.1749512 20554.18 0.1792613 1.102483 3.663 0.0015
## face - house 17_guess - 17_high Right 0.9632012 0.2094179 3906.37 0.4104213 1.515981 4.599 <.0001
## face - house 17_guess - 200_high Right 2.8676323 0.1862807 29077.11 2.3761421 3.359122 15.394 <.0001
## face - house 17_low - 17_high Right 0.3223288 0.1850160 5675.08 -0.1659622 0.810620 1.742 0.4892
## face - house 17_low - 200_high Right 2.2267599 0.1657007 37310.22 1.7895756 2.663944 13.438 <.0001
## face - house 17_high - 200_high Right 1.9044311 0.1902513 11186.98 1.4024099 2.406452 10.010 <.0001
##
## Degrees-of-freedom method: satterthwaite
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 6 estimates
## P value adjustment: bonferroni method for 6 tests
summary(contra_SC_N170_spec[7:12], delta = delta_null, level = .9)
summary(contra_SC_N170_spec[7:12], delta = delta_null, side = ">")
summary(contra_SC_N170_spec[7:12], delta = delta_null, side = "<")
Comparisons of face-specific N170 amplitude component lateralization for different durations and confidence:
contra_SC_N170_spec_hemi <- contrast(emmeans(lmm.opt.N170.SC.E205, ~ Category + DuraConf + Hemisphere), #
interaction = "pairwise", adjust = "Bonferroni", infer = TRUE)
contra_SC_N170_spec_hemi
## Category_pairwise DuraConf_pairwise Hemisphere_pairwise estimate SE df lower.CL upper.CL t.ratio p.value
## face - house 17_guess - 17_low Left - Right -0.2732653 0.2464975 17365.60 -0.9236637 0.3771332 -1.109 1.0000
## face - house 17_guess - 17_high Left - Right -0.2145503 0.2931618 2919.39 -0.9885142 0.5594137 -0.732 1.0000
## face - house 17_guess - 200_high Left - Right -0.4831245 0.2629621 27102.35 -1.1769371 0.2106881 -1.837 0.3971
## face - house 17_low - 17_high Left - Right 0.0587150 0.2597280 4296.86 -0.6268319 0.7442619 0.226 1.0000
## face - house 17_low - 200_high Left - Right -0.2098592 0.2341167 36563.69 -0.8275529 0.4078344 -0.896 1.0000
## face - house 17_high - 200_high Left - Right -0.2685742 0.2674890 8591.17 -0.9744424 0.4372940 -1.004 1.0000
##
## Degrees-of-freedom method: satterthwaite
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 6 estimates
## P value adjustment: bonferroni method for 6 tests
summary(contra_SC_N170_spec_hemi, delta = delta_null, level = .9)
summary(contra_SC_N170_spec_hemi, delta = delta_null, side = ">")
summary(contra_SC_N170_spec_hemi, delta = delta_null, side = "<")
Other results
Comparisons of N170 amplitudes for different durations and confidence (without multiple comparison corrections):
summary(contra_SC_N170[1:24], adjust = "none")
summary(contra_SC_N170[1:24], adjust = "none", delta = delta_null, level = .9)
Linear mixed model with high subjective confidences
tmpBlock.count <- {
clean_beha_205 %>%
filter(Type == "intact") %>%
group_by(SubjCode, Category, Duration, blockID) %>%
summarize(Count_sum = n())
}
avg.resp.long.E205.Block <- {
clean_beha_205 %>%
filter(Type == "intact") %>%
group_by(SubjCode, Category, Duration, blockID, Resp) %>%
summarize(Count = n(), .groups = "drop") %>%
right_join(., tidyr::expand(., SubjCode, blockID, Category, Duration)) %>%
mutate(Count = if_else(is.na(Count), 0L, Count)) %>%
right_join(tmpBlock.count, by = c("SubjCode", "blockID", "Category", "Duration")) %>%
mutate(RespRate = Count / Count_sum) %>%
ungroup()
}
# descriptive statistics of accuracy for plotting
desc.resp.E205.Block <- {
avg.resp.long.E205.Block %>%
group_by(blockID, Category, Duration, Resp) %>%
summarize(Mean = mean(RespRate)
, N = n()
, SD = sd(RespRate)
, SE = SD/sqrt(N)
, SE.lo = Mean - SE, SE.hi = Mean + SE
, CI = SE * qt(0.975,N)
, Median = median(RespRate)
, Lower = Mean - SD, Upper = Mean + SD # for rainclound plot
) %>%
ungroup()
}
avg.resp.wide.E205.Block <- {
avg.resp.long.E205.Block %>%
mutate(Conditions = paste(blockID, Category, Duration, Resp, sep = "_")) %>%
select(SubjCode, Conditions, RespRate) %>%
spread(Conditions, RespRate)
}
block.high <- {
avg.resp.long.E205.Block %>%
filter(Category == "face", Duration == 17, Resp == 1, SubjCode != "P513")
}
t.test(filter(block.high, blockID == 1)$RespRate, filter(block.high, blockID == 2)$RespRate, paired = T)
##
## Paired t-test
##
## data: filter(block.high, blockID == 1)$RespRate and filter(block.high, blockID == 2)$RespRate
## t = -0.64323, df = 26, p-value = 0.5257
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.10424318 0.05455182
## sample estimates:
## mean of the differences
## -0.02484568
knitr::kable(desc.resp.E205.Block, digits = 4)
| 1 |
face |
17 |
1 |
0.4628 |
28 |
0.3419 |
0.0646 |
0.3982 |
0.5274 |
0.1324 |
0.3812 |
0.1209 |
0.8047 |
| 1 |
face |
17 |
2 |
0.3192 |
28 |
0.2240 |
0.0423 |
0.2769 |
0.3615 |
0.0867 |
0.3354 |
0.0952 |
0.5432 |
| 1 |
face |
17 |
3 |
0.1737 |
28 |
0.1856 |
0.0351 |
0.1386 |
0.2087 |
0.0719 |
0.0958 |
-0.0120 |
0.3593 |
| 1 |
face |
17 |
4 |
0.0472 |
25 |
0.0638 |
0.0128 |
0.0344 |
0.0599 |
0.0263 |
0.0208 |
-0.0166 |
0.1109 |
| 1 |
face |
17 |
5 |
0.0156 |
4 |
0.0202 |
0.0101 |
0.0055 |
0.0257 |
0.0281 |
0.0063 |
-0.0046 |
0.0359 |
| 1 |
house |
17 |
1 |
0.0120 |
17 |
0.0107 |
0.0026 |
0.0094 |
0.0146 |
0.0055 |
0.0083 |
0.0013 |
0.0227 |
| 1 |
house |
17 |
2 |
0.0658 |
28 |
0.0806 |
0.0152 |
0.0505 |
0.0810 |
0.0312 |
0.0417 |
-0.0148 |
0.1464 |
| 1 |
house |
17 |
3 |
0.4188 |
28 |
0.2770 |
0.0524 |
0.3664 |
0.4711 |
0.1072 |
0.3542 |
0.1417 |
0.6958 |
| 1 |
house |
17 |
4 |
0.3890 |
27 |
0.2334 |
0.0449 |
0.3441 |
0.4340 |
0.0922 |
0.4167 |
0.1556 |
0.6225 |
| 1 |
house |
17 |
5 |
0.1961 |
19 |
0.2386 |
0.0547 |
0.1413 |
0.2508 |
0.1146 |
0.1042 |
-0.0425 |
0.4347 |
| 2 |
face |
17 |
1 |
0.5046 |
27 |
0.3326 |
0.0640 |
0.4406 |
0.5686 |
0.1313 |
0.4750 |
0.1721 |
0.8372 |
| 2 |
face |
17 |
2 |
0.3786 |
24 |
0.2660 |
0.0543 |
0.3243 |
0.4329 |
0.1121 |
0.3875 |
0.1126 |
0.6446 |
| 2 |
face |
17 |
3 |
0.1300 |
25 |
0.1197 |
0.0239 |
0.1061 |
0.1539 |
0.0493 |
0.0750 |
0.0103 |
0.2497 |
| 2 |
face |
17 |
4 |
0.0863 |
21 |
0.0828 |
0.0181 |
0.0682 |
0.1044 |
0.0376 |
0.0500 |
0.0035 |
0.1691 |
| 2 |
face |
17 |
5 |
0.0225 |
10 |
0.0184 |
0.0058 |
0.0167 |
0.0283 |
0.0130 |
0.0125 |
0.0041 |
0.0409 |
| 2 |
face |
200 |
1 |
0.9768 |
28 |
0.0320 |
0.0061 |
0.9707 |
0.9828 |
0.0124 |
0.9875 |
0.9447 |
1.0088 |
| 2 |
face |
200 |
2 |
0.0400 |
10 |
0.0403 |
0.0127 |
0.0273 |
0.0527 |
0.0284 |
0.0188 |
-0.0003 |
0.0803 |
| 2 |
face |
200 |
3 |
0.0208 |
9 |
0.0125 |
0.0042 |
0.0167 |
0.0250 |
0.0094 |
0.0125 |
0.0083 |
0.0333 |
| 2 |
face |
200 |
4 |
0.0125 |
1 |
NA |
NA |
NA |
NA |
NA |
0.0125 |
NA |
NA |
| 2 |
face |
200 |
5 |
0.0125 |
4 |
0.0000 |
0.0000 |
0.0125 |
0.0125 |
0.0000 |
0.0125 |
0.0125 |
0.0125 |
| 2 |
house |
17 |
1 |
0.0156 |
4 |
0.0063 |
0.0031 |
0.0125 |
0.0187 |
0.0087 |
0.0125 |
0.0094 |
0.0219 |
| 2 |
house |
17 |
2 |
0.0354 |
18 |
0.0412 |
0.0097 |
0.0257 |
0.0451 |
0.0204 |
0.0250 |
-0.0058 |
0.0766 |
| 2 |
house |
17 |
3 |
0.3032 |
27 |
0.2482 |
0.0478 |
0.2555 |
0.3510 |
0.0980 |
0.3000 |
0.0550 |
0.5515 |
| 2 |
house |
17 |
4 |
0.4696 |
28 |
0.2377 |
0.0449 |
0.4247 |
0.5146 |
0.0920 |
0.4563 |
0.2319 |
0.7074 |
| 2 |
house |
17 |
5 |
0.2208 |
27 |
0.2831 |
0.0545 |
0.1664 |
0.2753 |
0.1118 |
0.0750 |
-0.0622 |
0.5039 |
| 2 |
house |
200 |
1 |
0.0187 |
6 |
0.0105 |
0.0043 |
0.0145 |
0.0230 |
0.0104 |
0.0125 |
0.0083 |
0.0292 |
| 2 |
house |
200 |
3 |
0.0150 |
10 |
0.0079 |
0.0025 |
0.0125 |
0.0175 |
0.0056 |
0.0125 |
0.0071 |
0.0229 |
| 2 |
house |
200 |
4 |
0.0331 |
20 |
0.0208 |
0.0046 |
0.0285 |
0.0378 |
0.0097 |
0.0250 |
0.0123 |
0.0539 |
| 2 |
house |
200 |
5 |
0.9670 |
28 |
0.0262 |
0.0049 |
0.9620 |
0.9719 |
0.0101 |
0.9750 |
0.9408 |
0.9931 |
rt.RainPlot.E205.Block <- {
ggplot(data = avg.resp.long.E205.Block, aes(y = RespRate, x = as.factor(Resp), fill = Category)) +
geom_flat_violin(position = position_nudge(x = .2, y = 0), alpha = .7) +
geom_point(aes(y = RespRate, color = Category), position = position_jitter(width = .15), size = .5, alpha = .8) +
geom_boxplot(aes(y = RespRate), width = .2, outlier.shape = NA, alpha = .7) +
geom_point(data = desc.resp.E205.Block, aes(y = Mean, x = Resp, color = Category), position = position_nudge(x = 0.3), size = 2.5) +
geom_errorbar(data = desc.resp.E205.Block, aes(y = Mean, ymin = Lower, ymax = Upper), position = position_nudge(x = 0.3), width = 0) +
facet_grid(Duration ~ blockID) +
# facet_grid(. ~ Type, switch = "x") + # create two panels to show data
# scale_colour_grey() + # start = .1, end = .6, color for the contour
# scale_fill_grey() + # start = .3, end = .6, color for the fill
# scale_color_brewer(palette = "Set1") + # palette = "RdBu"
# scale_fill_brewer(palette = "Set1") + # palette = "RdBu"
labs(title = "Reaction times for E205", x = "Stimulus Type", y = "Reaction times (ms)", fill = "Duration(ms)", color = "Duration(ms)") + # set the names for main, x and y axises and the legend
# coord_cartesian(ylim = c(0, 1.05)) + # set the limit for y axis
# scale_y_continuous(expand= c(0, 0)) + # remove the space between columns and x axis
theme_bw() + # the backgroud color
plot_theme +
NULL
}
Â
rt.RainPlot.E205.Block
bloc.labels <- c("1st Half", "2nd Half")
names(bloc.labels) <- c(1, 2)
acc.ColuPlot.E205.Block <- {
ggplot(data = desc.resp.E205.Block, aes(y = Mean, x = Resp, fill = Category)) + # set the data, varialbes for x and y axises, and the variable for dividing data
geom_col(position = "dodge", color = "black", alpha = .7) + # position of columns and countour of columns
facet_grid(Duration ~ blockID,
labeller = labeller(Duration = dura.labels,
blockID = bloc.labels)) +
geom_errorbar(mapping = aes(ymin = SE.lo, ymax = SE.hi), linetype = 1, # set the error bar
show.legend = FALSE, width = 0.25, alpha = .5,
position = position_dodge(width=0.9)) +
geom_hline(yintercept = c(1), linetype = 5, alpha = 0.5) + # add the line for 0.5 and 1 (y)
coord_cartesian(ylim = c(0,1.1)) + # set the limit for y axis
scale_y_continuous(expand= c(0, 0)) + # remove the space between columns and x axis
labs(title = "", x = "Responses", y = "Percentage", fill = "Category") + # set the names for main, x and y axises Subjective Responses for E205
scale_fill_grey() + # set the color for columns
# geom_text(label = c("***", "***", "", ""), color = c("red", "red", "red", "red"), size = 6, nudge_y = 0.05, nudge_x = rep(c(-0.225, 0.225), 2)) + # add starts to the significant columns
erp_theme
}
acc.ColuPlot.E205.Block
df.ratio.E205.Block <- {
clean_beha_205 %>%
filter(Type == "intact") %>%
mutate(isFace = if_else(Category == "face", 1, 0)) %>%
group_by(SubjCode, blockID, Duration, Resp) %>%
summarize(ratio = mean(isFace), Count = n(), .groups = "drop")
}
desc.ratio.E205.Block <- {
df.ratio.E205.Block %>%
group_by(blockID, Duration, Resp) %>%
summarize(Mean = mean(ratio)
, N = n()
, SD = sd(ratio)
, SE = SD/sqrt(N)
, SE.lo = Mean - SE, SE.hi = Mean + SE
, CI = SE * qt(0.975,N)
, Median = median(ratio)
, Lower = Mean - SD, Upper = Mean + SD # for rainclound plot
) %>%
ungroup()
}
ratio.ColuPlot.E205.Block <- {
ggplot(data = desc.ratio.E205.Block, aes(y = Mean, x = Resp)) + # set the data, varialbes for x and y axises, and the variable for dividing data
geom_col(position = "dodge", color = "black", alpha = .7) + # position of columns and countour of columns
facet_grid(Duration ~ blockID,
labeller = labeller(Duration = dura.labels,
blockID = bloc.labels)) +
geom_errorbar(mapping = aes(ymin = SE.lo, ymax = SE.hi), linetype = 1, # set the error bar
show.legend = FALSE, width = 0.25, alpha = .5,
position = position_dodge(width=0.9)) +
geom_hline(yintercept = c(0.5, 1), linetype = 5, alpha = 0.5) + # add the line for 0.5 and 1 (y)
coord_cartesian(ylim = c(0,1.1)) + # set the limit for y axis
scale_y_continuous(expand= c(0, 0)) + # remove the space between columns and x axis
labs(title = "", x = "Responses", y = "Percentage") + # set the names for main, x and y axises Subjective Responses for E205
scale_fill_grey() + # set the color for columns
# geom_text(label = c("***", "***", "", ""), color = c("red", "red", "red", "red"), size = 6, nudge_y = 0.05, nudge_x = rep(c(-0.225, 0.225), 2)) + # add starts to the significant columns
plot_theme
}
ratio.ColuPlot.E205.Block
df.erp.high.E2 <- {
df.erp.SC.E2 %>%
filter(Confidence == "high") %>%
mutate(DuraBlock = paste(Duration, Block_ST_table, sep = "_")) %>%
sdif_coding_E205_erp_conf()
}
Amplitudes of the P1
# only keep the data for amplitudes of the P1 (correct and incorrect erps)
df.erp.P1.high.E2 <- {
df.erp.high.E2 %>%
filter(Component == "P1")
}
if (saveCSV) {
output.erp.P1.high.E2 <- file.path(folder_nesi, "E205_erp_P1_high.RData")
save(df.erp.P1.high.E2, file = output.erp.P1.high.E2)
}
# the structure of this dataset
head(df.erp.P1.high.E2, 10)
The maximal model
lmm.max.P1.high.E2 <- lmer(MeanAmp ~ Category * Hemisphere * DuraBlock +
(1 + Cate_C + Hemi_C + Cate_Hemi | SubjCode) +
(1 + Hemi_C | Stimuli),
data = df.erp.P1.high.E2,
verbose = TRUE,
control = lmerControl(optimizer = "bobyqa", # nloptwrap Nelder_Mead
optCtrl = list(maxfun = 1e5)))
## start par. = 1 0 1 1 0 0 0 1 0 0 1 0 1 fn = 104824.6
## At return
## eval: 865 fn: 103991.60 par: 0.0607829 0.0132275 0.00000 0.368838 -0.0559286 0.0565163 0.0161717 0.107212 0.0588553 -0.00506727 0.178642 -0.0172575 0.00000
print(summary(lmm.max.P1.high.E2), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Category * Hemisphere * DuraBlock + (1 + Cate_C + Hemi_C + Cate_Hemi | SubjCode) + (1 + Hemi_C | Stimuli)
## Data: df.erp.P1.high.E2
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+05))
##
## REML criterion at convergence: 103991.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.3462 -0.5802 -0.0117 0.5811 6.3206
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Stimuli (Intercept) 0.059710 0.24436
## Hemi_C 0.002828 0.05318 1.00
## SubjCode (Intercept) 2.198631 1.48278
## Cate_C 0.236320 0.48613 -0.46
## Hemi_C 0.623365 0.78953 0.29 0.13
## Cate_Hemi 0.009455 0.09724 0.67 -0.49 -0.52
## Residual 16.161464 4.02013
## Number of obs: 18456, groups: Stimuli, 80; SubjCode, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.379e+00 2.848e-01 2.794e+01 4.843 4.29e-05 ***
## Category2-1 -3.119e-03 1.354e-01 4.647e+01 -0.023 0.982
## Hemisphere2-1 1.101e-01 1.695e-01 3.103e+01 0.650 0.521
## DuraBlock2-1 -3.110e-02 9.951e-02 1.778e+04 -0.313 0.755
## DuraBlock3-2 1.062e-01 9.877e-02 1.193e+04 1.075 0.282
## Category2-1:Hemisphere2-1 -1.317e-01 1.491e-01 1.563e+03 -0.883 0.377
## Category2-1:DuraBlock2-1 -2.797e-02 1.988e-01 1.724e+04 -0.141 0.888
## Category2-1:DuraBlock3-2 1.071e-01 1.947e-01 5.630e+03 0.550 0.582
## Hemisphere2-1:DuraBlock2-1 4.922e-02 1.975e-01 1.801e+04 0.249 0.803
## Hemisphere2-1:DuraBlock3-2 1.643e-02 1.924e-01 1.163e+04 0.085 0.932
## Category2-1:Hemisphere2-1:DuraBlock2-1 -1.170e-02 3.940e-01 1.818e+04 -0.030 0.976
## Category2-1:Hemisphere2-1:DuraBlock3-2 -1.721e-01 3.737e-01 1.591e+04 -0.461 0.645
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## optimizer (bobyqa) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
The zero-correlation-parameter model
lmm.zcp.P1.high.E2 <- update(lmm.max.P1.high.E2,
formula = MeanAmp ~ Category * Hemisphere * DuraBlock +
(1 + Cate_C + Hemi_C + Cate_Hemi || SubjCode) +
(1 + Hemi_C || Stimuli))
## start par. = 1 1 1 1 1 1 fn = 104824.6
## At return
## eval: 349 fn: 103999.59 par: 0.00000 0.0606616 0.00000 0.196283 0.126125 0.371879
print(summary(lmm.zcp.P1.high.E2), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Category + Hemisphere + DuraBlock + (1 + Cate_C + Hemi_C + Cate_Hemi || SubjCode) + (1 + Hemi_C || Stimuli) + Category:Hemisphere + Category:DuraBlock + Hemisphere:DuraBlock + Category:Hemisphere:DuraBlock
## Data: df.erp.P1.high.E2
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+05))
##
## REML criterion at convergence: 103999.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.3295 -0.5808 -0.0128 0.5814 6.3179
##
## Random effects:
## Groups Name Variance Std.Dev.
## Stimuli Hemi_C 0.00000 0.0000
## Stimuli.1 (Intercept) 0.05947 0.2439
## SubjCode Cate_Hemi 0.00000 0.0000
## SubjCode.1 Hemi_C 0.62267 0.7891
## SubjCode.2 Cate_C 0.25710 0.5070
## SubjCode.3 (Intercept) 2.23512 1.4950
## Residual 16.16202 4.0202
## Number of obs: 18456, groups: Stimuli, 80; SubjCode, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.382e+00 2.872e-01 2.798e+01 4.811 4.66e-05 ***
## Category2-1 6.224e-03 1.392e-01 4.769e+01 0.045 0.965
## Hemisphere2-1 1.018e-01 1.696e-01 3.174e+01 0.601 0.552
## DuraBlock2-1 -3.233e-02 9.964e-02 1.821e+04 -0.324 0.746
## DuraBlock3-2 1.035e-01 9.925e-02 1.579e+04 1.043 0.297
## Category2-1:Hemisphere2-1 -1.448e-01 1.482e-01 1.785e+04 -0.977 0.328
## Category2-1:DuraBlock2-1 -2.736e-02 1.991e-01 1.794e+04 -0.137 0.891
## Category2-1:DuraBlock3-2 9.423e-02 1.962e-01 9.184e+03 0.480 0.631
## Hemisphere2-1:DuraBlock2-1 4.850e-02 1.976e-01 1.831e+04 0.245 0.806
## Hemisphere2-1:DuraBlock3-2 2.796e-02 1.929e-01 1.636e+04 0.145 0.885
## Category2-1:Hemisphere2-1:DuraBlock2-1 -6.056e-03 3.941e-01 1.832e+04 -0.015 0.988
## Category2-1:Hemisphere2-1:DuraBlock3-2 -1.563e-01 3.742e-01 1.830e+04 -0.418 0.676
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## optimizer (bobyqa) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
The reduced model
summary(rePCA(lmm.zcp.P1.high.E2))
## $Stimuli
## Importance of components:
## [,1] [,2]
## Standard deviation 0.06066 0
## Proportion of Variance 1.00000 0
## Cumulative Proportion 1.00000 1
##
## $SubjCode
## Importance of components:
## [,1] [,2] [,3] [,4]
## Standard deviation 0.3719 0.1963 0.12613 0
## Proportion of Variance 0.7176 0.1999 0.08254 0
## Cumulative Proportion 0.7176 0.9175 1.00000 1
lmm.rdc.P1.high.E2 <- update(lmm.zcp.P1.high.E2,
formula = MeanAmp ~ Category * Hemisphere * DuraBlock +
(1 + Cate_C + Hemi_C || SubjCode) +
(1 | Stimuli))
## start par. = 1 1 1 1 fn = 104459.8
## At return
## eval: 150 fn: 103999.59 par: 0.0606616 0.196283 0.126125 0.371878
print(summary(lmm.rdc.P1.high.E2), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Category + Hemisphere + DuraBlock + (1 + Cate_C + Hemi_C || SubjCode) + (1 | Stimuli) + Category:Hemisphere + Category:DuraBlock + Hemisphere:DuraBlock + Category:Hemisphere:DuraBlock
## Data: df.erp.P1.high.E2
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+05))
##
## REML criterion at convergence: 103999.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.3295 -0.5808 -0.0128 0.5814 6.3179
##
## Random effects:
## Groups Name Variance Std.Dev.
## Stimuli (Intercept) 0.05947 0.2439
## SubjCode Hemi_C 0.62267 0.7891
## SubjCode.1 Cate_C 0.25710 0.5070
## SubjCode.2 (Intercept) 2.23510 1.4950
## Residual 16.16202 4.0202
## Number of obs: 18456, groups: Stimuli, 80; SubjCode, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.382e+00 2.872e-01 2.798e+01 4.811 4.66e-05 ***
## Category2-1 6.224e-03 1.392e-01 4.769e+01 0.045 0.965
## Hemisphere2-1 1.018e-01 1.696e-01 3.174e+01 0.601 0.552
## DuraBlock2-1 -3.233e-02 9.964e-02 1.821e+04 -0.324 0.746
## DuraBlock3-2 1.035e-01 9.925e-02 1.579e+04 1.043 0.297
## Category2-1:Hemisphere2-1 -1.448e-01 1.482e-01 1.785e+04 -0.977 0.328
## Category2-1:DuraBlock2-1 -2.736e-02 1.991e-01 1.794e+04 -0.137 0.891
## Category2-1:DuraBlock3-2 9.423e-02 1.962e-01 9.184e+03 0.480 0.631
## Hemisphere2-1:DuraBlock2-1 4.850e-02 1.976e-01 1.831e+04 0.245 0.806
## Hemisphere2-1:DuraBlock3-2 2.796e-02 1.929e-01 1.636e+04 0.145 0.885
## Category2-1:Hemisphere2-1:DuraBlock2-1 -6.056e-03 3.941e-01 1.832e+04 -0.015 0.988
## Category2-1:Hemisphere2-1:DuraBlock3-2 -1.563e-01 3.742e-01 1.830e+04 -0.418 0.676
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lmm.rdc.P1.high.E2, lmm.zcp.P1.high.E2, refit = FALSE)
The extended model
lmm.etd.P1.high.E2 <- update(lmm.rdc.P1.high.E2,
formula = MeanAmp ~ Category * Hemisphere * DuraBlock +
(1 + Cate_C + Hemi_C | SubjCode) +
(1 | Stimuli))
## start par. = 1 1 0 0 1 0 1 fn = 104459.8
## At return
## eval: 323 fn: 103992.90 par: 0.0606162 0.368822 -0.0560008 0.0540837 0.107121 0.0594280 0.180418
print(summary(lmm.etd.P1.high.E2), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Category + Hemisphere + DuraBlock + (1 + Cate_C + Hemi_C | SubjCode) + (1 | Stimuli) + Category:Hemisphere + Category:DuraBlock + Hemisphere:DuraBlock + Category:Hemisphere:DuraBlock
## Data: df.erp.P1.high.E2
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+05))
##
## REML criterion at convergence: 103992.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.3330 -0.5802 -0.0119 0.5800 6.3237
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Stimuli (Intercept) 0.05939 0.2437
## SubjCode (Intercept) 2.19869 1.4828
## Cate_C 0.23616 0.4860 -0.46
## Hemi_C 0.63049 0.7940 0.27 0.14
## Residual 16.16335 4.0204
## Number of obs: 18456, groups: Stimuli, 80; SubjCode, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.379e+00 2.848e-01 2.794e+01 4.842 4.29e-05 ***
## Category2-1 -3.143e-03 1.353e-01 4.641e+01 -0.023 0.982
## Hemisphere2-1 1.061e-01 1.703e-01 3.162e+01 0.623 0.538
## DuraBlock2-1 -3.146e-02 9.952e-02 1.778e+04 -0.316 0.752
## DuraBlock3-2 1.067e-01 9.878e-02 1.192e+04 1.081 0.280
## Category2-1:Hemisphere2-1 -1.444e-01 1.481e-01 1.772e+04 -0.975 0.330
## Category2-1:DuraBlock2-1 -2.871e-02 1.988e-01 1.724e+04 -0.144 0.885
## Category2-1:DuraBlock3-2 1.074e-01 1.947e-01 5.619e+03 0.552 0.581
## Hemisphere2-1:DuraBlock2-1 4.903e-02 1.976e-01 1.830e+04 0.248 0.804
## Hemisphere2-1:DuraBlock3-2 2.236e-02 1.928e-01 1.603e+04 0.116 0.908
## Category2-1:Hemisphere2-1:DuraBlock2-1 -9.847e-03 3.941e-01 1.832e+04 -0.025 0.980
## Category2-1:Hemisphere2-1:DuraBlock3-2 -1.557e-01 3.742e-01 1.829e+04 -0.416 0.677
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lmm.etd.P1.high.E2, lmm.rdc.P1.high.E2, refit = FALSE)
The optimal model
lmm.opt.P1.high.E205 <- lmm.etd.P1.high.E2
summary(lmm.opt.P1.high.E205)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Category + Hemisphere + DuraBlock + (1 + Cate_C + Hemi_C | SubjCode) + (1 | Stimuli) + Category:Hemisphere + Category:DuraBlock + Hemisphere:DuraBlock + Category:Hemisphere:DuraBlock
## Data: df.erp.P1.high.E2
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+05))
##
## REML criterion at convergence: 103992.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.3330 -0.5802 -0.0119 0.5800 6.3237
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Stimuli (Intercept) 0.05939 0.2437
## SubjCode (Intercept) 2.19869 1.4828
## Cate_C 0.23616 0.4860 -0.46
## Hemi_C 0.63049 0.7940 0.27 0.14
## Residual 16.16335 4.0204
## Number of obs: 18456, groups: Stimuli, 80; SubjCode, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.379e+00 2.848e-01 2.794e+01 4.842 4.29e-05 ***
## Category2-1 -3.143e-03 1.353e-01 4.641e+01 -0.023 0.982
## Hemisphere2-1 1.061e-01 1.703e-01 3.162e+01 0.623 0.538
## DuraBlock2-1 -3.146e-02 9.952e-02 1.778e+04 -0.316 0.752
## DuraBlock3-2 1.067e-01 9.878e-02 1.192e+04 1.081 0.280
## Category2-1:Hemisphere2-1 -1.444e-01 1.481e-01 1.772e+04 -0.975 0.330
## Category2-1:DuraBlock2-1 -2.871e-02 1.988e-01 1.724e+04 -0.144 0.885
## Category2-1:DuraBlock3-2 1.074e-01 1.947e-01 5.619e+03 0.552 0.581
## Hemisphere2-1:DuraBlock2-1 4.903e-02 1.976e-01 1.830e+04 0.248 0.804
## Hemisphere2-1:DuraBlock3-2 2.236e-02 1.928e-01 1.603e+04 0.116 0.908
## Category2-1:Hemisphere2-1:DuraBlock2-1 -9.847e-03 3.941e-01 1.832e+04 -0.025 0.980
## Category2-1:Hemisphere2-1:DuraBlock3-2 -1.557e-01 3.742e-01 1.829e+04 -0.416 0.677
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Ctg2-1 Hms2-1 DrB2-1 DrB3-2 Ct2-1:H2-1 C2-1:DB2 C2-1:DB3 H2-1:DB2 H2-1:DB3 C2-1:H2-1:DB2
## Category2-1 -0.267
## Hemisphr2-1 0.238 0.089
## DuraBlck2-1 0.036 0.027 -0.001
## DuraBlck3-2 -0.082 -0.193 -0.004 -0.660
## Ctg2-1:H2-1 0.001 0.005 0.184 -0.001 -0.004
## Ct2-1:DB2-1 0.007 0.158 -0.001 0.442 -0.252 -0.001
## Ct2-1:DB3-2 -0.047 -0.321 -0.005 -0.254 0.333 -0.004 -0.674
## Hm2-1:DB2-1 0.000 -0.001 0.131 0.001 0.001 0.070 0.001 0.001
## Hm2-1:DB3-2 -0.001 -0.007 -0.242 0.001 0.004 -0.279 0.001 0.007 -0.685
## C2-1:H2-1:DB2 0.000 0.000 0.033 0.001 0.000 0.316 0.000 0.000 0.439 -0.268
## C2-1:H2-1:DB3 -0.001 -0.002 -0.125 0.000 0.002 -0.468 0.000 0.002 -0.273 0.309 -0.714
anova(lmm.opt.P1.high.E205)
Diagnostic plots
qqplot_lmer(lmm.opt.P1.high.E205)
Estimated marginal means
emm.P1.high.E205 <- emmeans(lmm.opt.P1.high.E205, ~ DuraBlock + Category + Hemisphere)
emm.P1.high.E205 %>%
as_tibble() %>%
arrange(DuraBlock, Category, Hemisphere)
Plots
# emmip(emm.P1.high.E205, Category ~ DuraBlock | Hemisphere, CIs = TRUE)
durablock.labels = c("17ms\nhalf 1", "17ms\nhalf 2", "200ms\nhalf 2")
p1.high.LinePlot = {
ggplot(data = data.frame(emm.P1.high.E205), aes(y = emmean, x = DuraBlock, color = Category, group = Category)) + # set the data, varialbes for x and y axises, and the variable for dividing data
geom_point(size = 2) +
geom_line(aes(linetype = Category), size = .85) + # position = "dodge", alpha = .7
scale_linetype_manual(values=c("solid", "dashed")) +
geom_errorbar(mapping = aes(ymin = lower.CL, ymax = upper.CL), linetype = 1, # set the error bar
show.legend = FALSE, width = 0.25, alpha = .5) + # , position = position_dodge(width=0.9)
facet_grid(. ~ Hemisphere, # switch = "x",
labeller = labeller(Hemisphere = hemiLabel)) +
scale_y_reverse(limits = ylimit_p1) + # set the limit for y axis
scale_x_discrete(labels = durablock.labels) +
# scale_y_continuous(expand= c(0, 0)) + # remove the space between columns and x axis
labs(x = "Duration and Halves", y = expression(paste("Amplitudes (", mu, "V)")), fill = "Category") + # set the names for main, x and y axises title = "",
# scale_x_discrete(labels = c("left", "right")) +
# geom_text(label = c("", "", "", "", "***", "", "", ""), color = "red", size = 6, nudge_y = 0.5, nudge_x = 0.5) + # add starts to the significant columns
theme_bw() +
erp_theme
}
# p1.high.LinePlot.slide <- {
# p1.high.LinePlot +
# theme(legend.position="right",
# axis.line = element_line(color="black", size = 3))
# }
# ggsave('P1_high.png', p1.high.LinePlot.slide, width = 7.5, height = 4)
p1.high.LinePlot
duration_color <- c(brewer.pal(9, 'Blues')[c(6, 8)], brewer.pal(12, "Paired")[6])
# category_color <- brewer.pal(8, 'Dark2')[c(6,1)]
p1_high_LinePlot = {
data.frame(emm.P1.high.E205) %>%
mutate(DuraBlock = fct_recode(DuraBlock,
`33_1` = "17_1", `33_2` = "17_2", `216_2` = "200_2")) %>%
ggplot(aes(y = emmean, x = DuraBlock, group = Category)) + # set the data, varialbes for x and y axises, and the variable for dividing data
geom_point(aes(fill = DuraBlock), size=4, shape=21, stroke=0) +
scale_fill_manual(values = duration_color, guide = "none") +
geom_line(aes(linetype = Category)) + # position = "dodge", alpha = .7
scale_linetype_manual(values=c("solid", "dashed")) +
geom_errorbar(mapping = aes(ymin = lower.CL, ymax = upper.CL), linetype = 1, # set the error bar
show.legend = FALSE, width = 0.25, alpha = .5) + # , position = position_dodge(width=0.9)
facet_grid(. ~ Hemisphere, # switch = "x",
labeller = labeller(Hemisphere = hemiLabel)) +
scale_y_reverse(limits = ylimit_p1) + # set the limit for y axis
scale_x_discrete(labels = c("33 ms\nhalf 1", "33 ms\nhalf 2", "216 ms\nhalf 2")) +
# scale_y_continuous(expand= c(0, 0)) + # remove the space between columns and x axis
labs(x = "Duration and Halves", y = expression(paste("P1 Amplitudes (", mu, "V)")), fill = "Category") + # set the names for main, x and y axises title = "",
theme_bw() +
erp_theme +
theme(
axis.title.x = element_text(margin = margin(t = 10, r = 0, b = 0, l = 0)),
legend.margin = margin(t = -10),
legend.title=element_text(size=12),
legend.text=element_text(size=12),
)
}
ggsave('p1_high_pub.pdf', p1_high_LinePlot, width = 7, height = 4.5)
Contrast
Main results
Comparisons of P1 amplitudes for (intact) faces with different durations and high confidence (left hemisphere):
contra_high_P1 <- contrast(emmeans(lmm.opt.P1.high.E205, ~ DuraBlock | Hemisphere + Category),
"pairwise", adjust = "none", infer = TRUE)
contra_high_P1[1:3]
## contrast Hemisphere Category estimate SE df lower.CL upper.CL t.ratio p.value
## 17_1 - 17_2 Left face 0.04407997 0.1482603 18312.03 -0.2465240 0.3346840 0.297 0.7662
## 17_1 - 200_2 Left face 0.04115602 0.1227113 15818.16 -0.1993721 0.2816842 0.335 0.7373
## 17_2 - 200_2 Left face -0.00292394 0.1591442 17711.71 -0.3148623 0.3090144 -0.018 0.9853
##
## Degrees-of-freedom method: satterthwaite
## Confidence level used: 0.95
summary(contra_high_P1[1:3], delta = delta_null, level = .9)
summary(contra_high_P1[1:3], delta = delta_null, side = ">")
summary(contra_high_P1[1:3], delta = delta_null, side = "<")
Comparisons of P1 amplitudes for (intact) faces with different durations and high confidence (right hemisphere):
contra_high_P1[4:6]
## contrast Hemisphere Category estimate SE df lower.CL upper.CL t.ratio p.value
## 17_1 - 17_2 Right face -0.00987683 0.1481876 18277.67 -0.3003384 0.2805847 -0.067 0.9469
## 17_1 - 200_2 Right face -0.11299380 0.1223560 14516.06 -0.3528272 0.1268396 -0.923 0.3558
## 17_2 - 200_2 Right face -0.10311696 0.1589035 17271.69 -0.4145839 0.2083500 -0.649 0.5164
##
## Degrees-of-freedom method: satterthwaite
## Confidence level used: 0.95
summary(contra_high_P1[4:6], delta = delta_null, level = .9)
summary(contra_high_P1[4:6], delta = delta_null, side = ">")
summary(contra_high_P1[4:6], delta = delta_null, side = "<")
Comparisons of face-specific P1 amplitude components for (intact) faces with different durations and high confidence (left hemisphere):
contra_high_P1_spec <- contrast(emmeans(lmm.opt.P1.high.E205, ~ Category + DuraBlock | Hemisphere), #
interaction = "pairwise", adjust = "Bonferroni", infer = TRUE)
contra_high_P1_spec[1:3]
## Category_pairwise DuraBlock_pairwise Hemisphere estimate SE df lower.CL upper.CL t.ratio p.value
## face - house 17_1 - 17_2 Left -0.02378827 0.2798066 18069.83 -0.6937019 0.6461254 -0.085 1.0000
## face - house 17_1 - 200_2 Left 0.16146629 0.2148371 3745.20 -0.3530805 0.6760131 0.752 1.0000
## face - house 17_2 - 200_2 Left 0.18525456 0.2697033 10769.69 -0.4605105 0.8310196 0.687 1.0000
##
## Degrees-of-freedom method: satterthwaite
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 3 estimates
## P value adjustment: bonferroni method for 3 tests
summary(contra_high_P1_spec[1:3], delta = delta_null, level = .9)
summary(contra_high_P1_spec[1:3], delta = delta_null, side = ">")
summary(contra_high_P1_spec[1:3], delta = delta_null, side = "<")
Comparisons of face-specific P1 amplitude components for (intact) faces with different durations and high confidence (right hemisphere):
contra_high_P1_spec[4:6]
## Category_pairwise DuraBlock_pairwise Hemisphere estimate SE df lower.CL upper.CL t.ratio p.value
## face - house 17_1 - 17_2 Right -0.03363520 0.2799462 18184.19 -0.7038828 0.6366124 -0.120 1.0000
## face - house 17_1 - 200_2 Right -0.00405177 0.2160355 4863.00 -0.5214153 0.5133118 -0.019 1.0000
## face - house 17_2 - 200_2 Right 0.02958343 0.2703208 12065.72 -0.6176494 0.6768162 0.109 1.0000
##
## Degrees-of-freedom method: satterthwaite
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 3 estimates
## P value adjustment: bonferroni method for 3 tests
summary(contra_high_P1_spec[4:6], delta = delta_null, level = .9)
summary(contra_high_P1_spec[4:6], delta = delta_null, side = ">")
summary(contra_high_P1_spec[4:6], delta = delta_null, side = "<")
Other results
Comparisons of P1 amplitudes for different durations and high confidence (without multiple comparison corrections):
summary(contra_high_P1[1:12], adjust = "none")
summary(contra_high_P1[1:12], adjust = "none", delta = delta_null, level = .9)
Amplitudes of the N170
# only keep the data for amplitudes of the N170 (correct and incorrect erps)
df.erp.N170.high.E2 <- {
df.erp.high.E2 %>%
filter(Component == "N170")
}
if (saveCSV) {
output.erp.N170.high.E2 <- file.path(folder_nesi, "E205_erp_N170_high.RData")
save(df.erp.N170.high.E2, file = output.erp.N170.high.E2)
}
# the structure of this dataset
head(df.erp.N170.high.E2, 10)
The maximal model
lmm.max.N170.high.E2 <- lmer(MeanAmp ~ Category * Hemisphere * DuraBlock +
(1 + Cate_C + Hemi_C + Cate_Hemi | SubjCode) +
(1 + Hemi_C | Stimuli),
data = df.erp.N170.high.E2,
verbose = TRUE,
control = lmerControl(optimizer = "bobyqa", # nloptwrap Nelder_Mead
optCtrl = list(maxfun = 1e5)))
## start par. = 1 0 1 1 0 0 0 1 0 0 1 0 1 fn = 105474.1
## At return
## eval: 1238 fn: 104763.30 par: 0.0573010 0.0277339 0.00000 0.454977 -0.134782 -0.0159491 -0.0974498 0.235465 -0.0364216 -0.107640 0.347868 -0.233303 0.243926
print(summary(lmm.max.N170.high.E2), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Category * Hemisphere * DuraBlock + (1 + Cate_C + Hemi_C + Cate_Hemi | SubjCode) + (1 + Hemi_C | Stimuli)
## Data: df.erp.N170.high.E2
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+05))
##
## REML criterion at convergence: 104763.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.6334 -0.5935 -0.0012 0.5937 6.1298
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Stimuli (Intercept) 0.05502 0.2346
## Hemi_C 0.01289 0.1135 1.00
## SubjCode (Intercept) 3.46891 1.8625
## Cate_C 1.23353 1.1106 -0.50
## Hemi_C 2.05437 1.4333 -0.05 -0.07
## Cate_Hemi 2.26250 1.5042 -0.27 -0.12 -0.59
## Residual 16.75767 4.0936
## Number of obs: 18456, groups: Stimuli, 80; SubjCode, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -2.319e+00 3.558e-01 2.775e+01 -6.517 4.80e-07 ***
## Category2-1 1.625e+00 2.339e-01 3.373e+01 6.949 5.42e-08 ***
## Hemisphere2-1 -3.835e-01 2.846e-01 2.969e+01 -1.348 0.188
## DuraBlock2-1 -8.716e-01 1.016e-01 1.831e+04 -8.581 < 2e-16 ***
## DuraBlock3-2 1.019e+00 1.014e-01 1.710e+04 10.052 < 2e-16 ***
## Category2-1:Hemisphere2-1 4.113e-01 3.317e-01 3.519e+01 1.240 0.223
## Category2-1:DuraBlock2-1 -1.442e-01 2.030e-01 1.827e+04 -0.710 0.478
## Category2-1:DuraBlock3-2 1.909e+00 2.018e-01 1.457e+04 9.460 < 2e-16 ***
## Hemisphere2-1:DuraBlock2-1 -1.955e-01 2.024e-01 1.765e+04 -0.966 0.334
## Hemisphere2-1:DuraBlock3-2 4.096e-02 1.997e-01 1.025e+04 0.205 0.837
## Category2-1:Hemisphere2-1:DuraBlock2-1 4.030e-01 4.045e-01 1.739e+04 0.996 0.319
## Category2-1:Hemisphere2-1:DuraBlock3-2 3.516e-03 3.962e-01 6.953e+03 0.009 0.993
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## optimizer (bobyqa) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
The zero-correlation-parameter model
lmm.zcp.N170.high.E2 <- update(lmm.max.N170.high.E2,
formula = MeanAmp ~ Category * Hemisphere * DuraBlock +
(1 + Cate_C + Hemi_C + Cate_Hemi || SubjCode) +
(1 + Hemi_C || Stimuli))
## start par. = 1 1 1 1 1 1 fn = 105474.1
## At return
## eval: 474 fn: 104787.06 par: 8.11255e-08 0.0562735 0.374364 0.356541 0.272915 0.456968
print(summary(lmm.zcp.N170.high.E2), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Category + Hemisphere + DuraBlock + (1 + Cate_C + Hemi_C + Cate_Hemi || SubjCode) + (1 + Hemi_C || Stimuli) + Category:Hemisphere + Category:DuraBlock + Hemisphere:DuraBlock + Category:Hemisphere:DuraBlock
## Data: df.erp.N170.high.E2
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+05))
##
## REML criterion at convergence: 104787.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.6562 -0.5916 0.0000 0.5928 6.1351
##
## Random effects:
## Groups Name Variance Std.Dev.
## Stimuli Hemi_C 0.00000 0.0000
## Stimuli.1 (Intercept) 0.05308 0.2304
## SubjCode Cate_Hemi 2.34918 1.5327
## SubjCode.1 Hemi_C 2.13082 1.4597
## SubjCode.2 Cate_C 1.24849 1.1174
## SubjCode.3 (Intercept) 3.50025 1.8709
## Residual 16.76210 4.0942
## Number of obs: 18456, groups: Stimuli, 80; SubjCode, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -2.316e+00 3.574e-01 2.775e+01 -6.480 5.30e-07 ***
## Category2-1 1.639e+00 2.352e-01 3.368e+01 6.968 5.16e-08 ***
## Hemisphere2-1 -3.760e-01 2.898e-01 2.981e+01 -1.297 0.205
## DuraBlock2-1 -8.721e-01 1.017e-01 1.836e+04 -8.579 < 2e-16 ***
## DuraBlock3-2 1.015e+00 1.017e-01 1.820e+04 9.981 < 2e-16 ***
## Category2-1:Hemisphere2-1 4.367e-01 3.389e-01 3.544e+01 1.289 0.206
## Category2-1:DuraBlock2-1 -1.426e-01 2.032e-01 1.835e+04 -0.702 0.483
## Category2-1:DuraBlock3-2 1.888e+00 2.028e-01 1.705e+04 9.314 < 2e-16 ***
## Hemisphere2-1:DuraBlock2-1 -1.984e-01 2.029e-01 1.829e+04 -0.978 0.328
## Hemisphere2-1:DuraBlock3-2 3.172e-02 2.019e-01 1.734e+04 0.157 0.875
## Category2-1:Hemisphere2-1:DuraBlock2-1 4.145e-01 4.057e-01 1.823e+04 1.022 0.307
## Category2-1:Hemisphere2-1:DuraBlock3-2 -4.034e-02 4.019e-01 1.438e+04 -0.100 0.920
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## optimizer (bobyqa) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
The reduced model
summary(rePCA(lmm.zcp.N170.high.E2))
## $Stimuli
## Importance of components:
## [,1] [,2]
## Standard deviation 0.05627 0
## Proportion of Variance 1.00000 0
## Cumulative Proportion 1.00000 1
##
## $SubjCode
## Importance of components:
## [,1] [,2] [,3] [,4]
## Standard deviation 0.4570 0.3744 0.3565 0.2729
## Proportion of Variance 0.3793 0.2545 0.2309 0.1353
## Cumulative Proportion 0.3793 0.6338 0.8647 1.0000
lmm.rdc.N170.high.E2 <- update(lmm.zcp.N170.high.E2,
formula = MeanAmp ~ Category * Hemisphere * DuraBlock +
(1 + Cate_C + Hemi_C + Cate_Hemi || SubjCode) +
(1 | Stimuli))
## start par. = 1 1 1 1 1 fn = 105207.8
## At return
## eval: 259 fn: 104787.06 par: 0.0562734 0.374363 0.356541 0.272917 0.456967
print(summary(lmm.rdc.N170.high.E2), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Category + Hemisphere + DuraBlock + (1 + Cate_C + Hemi_C + Cate_Hemi || SubjCode) + (1 | Stimuli) + Category:Hemisphere + Category:DuraBlock + Hemisphere:DuraBlock + Category:Hemisphere:DuraBlock
## Data: df.erp.N170.high.E2
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+05))
##
## REML criterion at convergence: 104787.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.6562 -0.5916 0.0000 0.5928 6.1351
##
## Random effects:
## Groups Name Variance Std.Dev.
## Stimuli (Intercept) 0.05308 0.2304
## SubjCode Cate_Hemi 2.34917 1.5327
## SubjCode.1 Hemi_C 2.13082 1.4597
## SubjCode.2 Cate_C 1.24850 1.1174
## SubjCode.3 (Intercept) 3.50024 1.8709
## Residual 16.76210 4.0942
## Number of obs: 18456, groups: Stimuli, 80; SubjCode, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -2.316e+00 3.574e-01 2.775e+01 -6.480 5.30e-07 ***
## Category2-1 1.639e+00 2.352e-01 3.368e+01 6.968 5.16e-08 ***
## Hemisphere2-1 -3.760e-01 2.898e-01 2.981e+01 -1.297 0.205
## DuraBlock2-1 -8.721e-01 1.017e-01 1.836e+04 -8.579 < 2e-16 ***
## DuraBlock3-2 1.015e+00 1.017e-01 1.820e+04 9.981 < 2e-16 ***
## Category2-1:Hemisphere2-1 4.367e-01 3.389e-01 3.544e+01 1.289 0.206
## Category2-1:DuraBlock2-1 -1.426e-01 2.032e-01 1.835e+04 -0.702 0.483
## Category2-1:DuraBlock3-2 1.888e+00 2.028e-01 1.705e+04 9.314 < 2e-16 ***
## Hemisphere2-1:DuraBlock2-1 -1.984e-01 2.029e-01 1.829e+04 -0.978 0.328
## Hemisphere2-1:DuraBlock3-2 3.172e-02 2.019e-01 1.734e+04 0.157 0.875
## Category2-1:Hemisphere2-1:DuraBlock2-1 4.145e-01 4.057e-01 1.823e+04 1.022 0.307
## Category2-1:Hemisphere2-1:DuraBlock3-2 -4.034e-02 4.019e-01 1.438e+04 -0.100 0.920
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lmm.rdc.N170.high.E2, lmm.zcp.N170.high.E2, refit = FALSE)
The extended model
lmm.etd.N170.high.E2 <- update(lmm.rdc.N170.high.E2,
formula = MeanAmp ~ Category * Hemisphere * DuraBlock +
(1 + Cate_C + Hemi_C + Cate_Hemi | SubjCode) +
(1 | Stimuli))
## start par. = 1 1 0 0 0 1 0 0 1 0 1 fn = 105207.8
## At return
## eval: 903 fn: 104765.64 par: 0.0562480 0.454957 -0.134602 -0.0160593 -0.0976551 0.235449 -0.0361790 -0.107527 0.347685 -0.233521 0.243953
print(summary(lmm.etd.N170.high.E2), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Category + Hemisphere + DuraBlock + (1 + Cate_C + Hemi_C + Cate_Hemi | SubjCode) + (1 | Stimuli) + Category:Hemisphere + Category:DuraBlock + Hemisphere:DuraBlock + Category:Hemisphere:DuraBlock
## Data: df.erp.N170.high.E2
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+05))
##
## REML criterion at convergence: 104765.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.6367 -0.5923 -0.0010 0.5935 6.1423
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Stimuli (Intercept) 0.05303 0.2303
## SubjCode (Intercept) 3.46964 1.8627
## Cate_C 1.23296 1.1104 -0.50
## Hemi_C 2.05262 1.4327 -0.05 -0.07
## Cate_Hemi 2.26537 1.5051 -0.27 -0.12 -0.59
## Residual 16.76265 4.0942
## Number of obs: 18456, groups: Stimuli, 80; SubjCode, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -2.319e+00 3.558e-01 2.774e+01 -6.517 4.82e-07 ***
## Category2-1 1.626e+00 2.336e-01 3.361e+01 6.960 5.34e-08 ***
## Hemisphere2-1 -3.885e-01 2.842e-01 2.957e+01 -1.367 0.182
## DuraBlock2-1 -8.726e-01 1.016e-01 1.832e+04 -8.590 < 2e-16 ***
## DuraBlock3-2 1.019e+00 1.014e-01 1.712e+04 10.054 < 2e-16 ***
## Category2-1:Hemisphere2-1 4.041e-01 3.309e-01 3.480e+01 1.221 0.230
## Category2-1:DuraBlock2-1 -1.453e-01 2.031e-01 1.828e+04 -0.716 0.474
## Category2-1:DuraBlock3-2 1.908e+00 2.018e-01 1.459e+04 9.454 < 2e-16 ***
## Hemisphere2-1:DuraBlock2-1 -1.971e-01 2.024e-01 1.765e+04 -0.974 0.330
## Hemisphere2-1:DuraBlock3-2 5.056e-02 1.997e-01 1.025e+04 0.253 0.800
## Category2-1:Hemisphere2-1:DuraBlock2-1 3.992e-01 4.045e-01 1.739e+04 0.987 0.324
## Category2-1:Hemisphere2-1:DuraBlock3-2 1.533e-02 3.962e-01 6.960e+03 0.039 0.969
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lmm.etd.N170.high.E2, lmm.rdc.N170.high.E2, refit = FALSE)
The optimal model
lmm.opt.N170.high.E205 <- lmm.etd.N170.high.E2
summary(lmm.opt.N170.high.E205)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Category + Hemisphere + DuraBlock + (1 + Cate_C + Hemi_C + Cate_Hemi | SubjCode) + (1 | Stimuli) + Category:Hemisphere + Category:DuraBlock + Hemisphere:DuraBlock + Category:Hemisphere:DuraBlock
## Data: df.erp.N170.high.E2
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+05))
##
## REML criterion at convergence: 104765.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.6367 -0.5923 -0.0010 0.5935 6.1423
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Stimuli (Intercept) 0.05303 0.2303
## SubjCode (Intercept) 3.46964 1.8627
## Cate_C 1.23296 1.1104 -0.50
## Hemi_C 2.05262 1.4327 -0.05 -0.07
## Cate_Hemi 2.26537 1.5051 -0.27 -0.12 -0.59
## Residual 16.76265 4.0942
## Number of obs: 18456, groups: Stimuli, 80; SubjCode, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -2.319e+00 3.558e-01 2.774e+01 -6.517 4.82e-07 ***
## Category2-1 1.626e+00 2.336e-01 3.361e+01 6.960 5.34e-08 ***
## Hemisphere2-1 -3.885e-01 2.842e-01 2.957e+01 -1.367 0.182
## DuraBlock2-1 -8.726e-01 1.016e-01 1.832e+04 -8.590 < 2e-16 ***
## DuraBlock3-2 1.019e+00 1.014e-01 1.712e+04 10.054 < 2e-16 ***
## Category2-1:Hemisphere2-1 4.041e-01 3.309e-01 3.480e+01 1.221 0.230
## Category2-1:DuraBlock2-1 -1.453e-01 2.031e-01 1.828e+04 -0.716 0.474
## Category2-1:DuraBlock3-2 1.908e+00 2.018e-01 1.459e+04 9.454 < 2e-16 ***
## Hemisphere2-1:DuraBlock2-1 -1.971e-01 2.024e-01 1.765e+04 -0.974 0.330
## Hemisphere2-1:DuraBlock3-2 5.056e-02 1.997e-01 1.025e+04 0.253 0.800
## Category2-1:Hemisphere2-1:DuraBlock2-1 3.992e-01 4.045e-01 1.739e+04 0.987 0.324
## Category2-1:Hemisphere2-1:DuraBlock3-2 1.533e-02 3.962e-01 6.960e+03 0.039 0.969
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Ctg2-1 Hms2-1 DrB2-1 DrB3-2 Ct2-1:H2-1 C2-1:DB2 C2-1:DB3 H2-1:DB2 H2-1:DB3 C2-1:H2-1:DB2
## Category2-1 -0.418
## Hemisphr2-1 -0.044 -0.061
## DuraBlck2-1 0.029 0.013 0.001
## DuraBlck3-2 -0.070 -0.126 0.004 -0.651
## Ctg2-1:H2-1 -0.228 -0.101 -0.415 0.002 0.009
## Ct2-1:DB2-1 0.004 0.089 0.001 0.444 -0.245 0.002
## Ct2-1:DB3-2 -0.041 -0.207 0.005 -0.246 0.345 0.011 -0.656
## Hm2-1:DB2-1 0.000 0.001 0.075 -0.002 -0.001 0.024 -0.002 -0.001
## Hm2-1:DB3-2 0.002 0.006 -0.162 -0.001 -0.009 -0.154 -0.002 -0.010 -0.667
## C2-1:H2-1:DB2 0.000 0.002 0.014 -0.002 -0.002 0.131 -0.003 -0.002 0.440 -0.254
## C2-1:H2-1:DB3 0.002 0.008 -0.090 -0.001 -0.010 -0.267 -0.002 -0.013 -0.255 0.325 -0.674
anova(lmm.opt.N170.high.E205)
Diagnostic plots
qqplot_lmer(lmm.opt.N170.high.E205)
Estimated marginal means
emm.N170.high.E205 <- emmeans(lmm.opt.N170.high.E205, ~ DuraBlock + Category + Hemisphere)
emm.N170.high.E205 %>%
as_tibble() %>%
arrange(DuraBlock, Category, Hemisphere)
Plots
# emmip(emm.N170.high.E205, Category ~ DuraBlock | Hemisphere, CIs = TRUE) +
# scale_y_continuous(trans = "reverse")
n170.high.LinePlot = {
ggplot(data = data.frame(emm.N170.high.E205), aes(y = emmean, x = DuraBlock, color = Category, group = Category)) + # set the data, varialbes for x and y axises, and the variable for dividing data
geom_point(size = 2) +
geom_line(aes(linetype = Category), size = .85) + # position = "dodge", alpha = .7
scale_linetype_manual(values=c("solid", "dashed")) +
geom_errorbar(mapping = aes(ymin = lower.CL, ymax = upper.CL), linetype = 1, # set the error bar
show.legend = FALSE, width = 0.25, alpha = .5) + # , position = position_dodge(width=0.9)
facet_grid(. ~ Hemisphere, # switch = "x",
labeller = labeller(Hemisphere = hemiLabel)) +
scale_y_reverse(limits = ylimit_n170, breaks = seq(1, -5, -1)) + # set the limit for y axis
scale_x_discrete(labels = durablock.labels) +
# scale_y_continuous(expand= c(0, 0)) + # remove the space between columns and x axis
labs(x = "Duration and Halves", y = expression(paste("Amplitudes (", mu, "V)")), fill = "Category") + # set the names for main, x and y axises title = "",
# scale_x_discrete(labels = c("left", "right")) +
# geom_text(label = c("", "", "", "", "***", "", "", ""), color = "red", size = 6, nudge_y = 0.5, nudge_x = 0.5) + # add starts to the significant columns
theme_bw() +
erp_theme
}
# n170.high.LinePlot.slide <- {
# n170.high.LinePlot +
# scale_y_reverse(limits = c(0.8, -4.8), breaks = seq(0, -4, -1)) + # set the limit for y axis
# theme(legend.position="right",
# axis.line = element_line(color="black", size = 3))
# }
# ggsave('N170_high.png', n170.high.LinePlot.slide, width = 7.5, height = 4)
n170.high.LinePlot
duration_color <- c(brewer.pal(9, 'Blues')[c(6, 8)], brewer.pal(12, "Paired")[6])
# category_color <- brewer.pal(8, 'Dark2')[c(6,1)]
n170_high_LinePlot = {
data.frame(emm.N170.high.E205) %>%
mutate(DuraBlock = fct_recode(DuraBlock,
`33_1` = "17_1", `33_2` = "17_2", `216_2` = "200_2")) %>%
ggplot(aes(y = emmean, x = DuraBlock, group = Category)) + # set the data, varialbes for x and y axises, and the variable for dividing data
geom_point(aes(fill = DuraBlock), size=4, shape=21, stroke=0) +
scale_fill_manual(values = duration_color, guide = "none") +
geom_line(aes(linetype = Category)) + # position = "dodge", alpha = .7
scale_linetype_manual(values=c("solid", "dashed")) +
geom_errorbar(mapping = aes(ymin = lower.CL, ymax = upper.CL), linetype = 1, # set the error bar
show.legend = FALSE, width = 0.25, alpha = .5) + # , position = position_dodge(width=0.9)
facet_grid(. ~ Hemisphere, # switch = "x",
labeller = labeller(Hemisphere = hemiLabel)) +
scale_y_reverse(limits = ylimit_n170, breaks = seq(1, -5, -1)) + # set the limit for y axis
scale_x_discrete(labels = c("33 ms\nhalf 1", "33 ms\nhalf 2", "216 ms\nhalf 2")) +
# scale_y_continuous(expand= c(0, 0)) + # remove the space between columns and x axis
labs(x = "Duration and Halves", y = expression(paste("N170 Amplitudes (", mu, "V)")), fill = "Category") + # set the names for main, x and y axises title = "",
theme_bw() +
erp_theme +
theme(
axis.title.x = element_text(margin = margin(t = 10, r = 0, b = 0, l = 0)),
legend.margin = margin(t = -10),
legend.title=element_text(size=12),
legend.text=element_text(size=12),
)
}
ggsave('n170_high_pub.pdf', n170_high_LinePlot, width = 7, height = 4.5)
Contrast
Main results
Comparisons of N170 amplitudes for (intact) faces with different durations and high confidence (left hemisphere):
contra_high_N170 <- contrast(emmeans(lmm.opt.N170.high.E205, ~ DuraBlock | Hemisphere + Category), #
"pairwise", adjust = "Bonferroni", infer = TRUE)
# left hemisphere
contra_high_N170[1:3]
## contrast Hemisphere Category estimate SE df lower.CL upper.CL t.ratio p.value
## 17_1 - 17_2 Left face 0.6015598 0.1514381 18306.67 0.2389867 0.9641330 3.972 0.0002
## 17_1 - 200_2 Left face 0.5576724 0.1272494 17686.68 0.2530109 0.8623339 4.383 <.0001
## 17_2 - 200_2 Left face -0.0438874 0.1636515 18181.18 -0.4357021 0.3479273 -0.268 1.0000
##
## Degrees-of-freedom method: satterthwaite
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 3 estimates
## P value adjustment: bonferroni method for 3 tests
summary(contra_high_N170[1:3], delta = delta_null, level = .9)
summary(contra_high_N170[1:3], delta = delta_null, side = ">")
summary(contra_high_N170[1:3], delta = delta_null, side = "<")
Comparisons of N170 amplitudes for (intact) faces with different durations and high confidence (right hemisphere):
# right hemisphere
contra_high_N170[4:6]
## contrast Hemisphere Category estimate SE df lower.CL upper.CL t.ratio p.value
## 17_1 - 17_2 Right face 0.9982659 0.1514036 18303.96 0.6357754 1.3607564 6.593 <.0001
## 17_1 - 200_2 Right face 0.9114827 0.1270620 17334.05 0.6072693 1.2156961 7.174 <.0001
## 17_2 - 200_2 Right face -0.0867832 0.1635189 18079.17 -0.4782807 0.3047143 -0.531 1.0000
##
## Degrees-of-freedom method: satterthwaite
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 3 estimates
## P value adjustment: bonferroni method for 3 tests
summary(contra_high_N170[4:6], delta = delta_null, level = .9)
summary(contra_high_N170[4:6], delta = delta_null, side = ">")
summary(contra_high_N170[4:6], delta = delta_null, side = "<")
Comparisons of N170 amplitude lateralization for (intact) faces with different durations and high confidence:
contra_high_N170_hemi <- contrast(emmeans(lmm.opt.N170.high.E205, ~ DuraBlock + Hemisphere | Category), #
interaction = "pairwise", adjust = "Bonferroni", infer = TRUE)
contra_high_N170_hemi[1:3]
## DuraBlock_pairwise Hemisphere_pairwise Category estimate SE df lower.CL upper.CL t.ratio p.value
## 17_1 - 17_2 Left - Right face -0.3967060 0.2140373 18283.98 -0.9091541 0.1157421 -1.853 0.1915
## 17_1 - 200_2 Left - Right face -0.3538102 0.1794494 16921.73 -0.7834513 0.0758308 -1.972 0.1460
## 17_2 - 200_2 Left - Right face 0.0428958 0.2310594 17949.51 -0.5103077 0.5960992 0.186 1.0000
##
## Degrees-of-freedom method: satterthwaite
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 3 estimates
## P value adjustment: bonferroni method for 3 tests
summary(contra_high_N170_hemi[1:3], delta = delta_null, level = .9)
summary(contra_high_N170_hemi[1:3], delta = delta_null, side = ">")
summary(contra_high_N170_hemi[1:3], delta = delta_null, side = "<")
Comparisons of face-specific N170 amplitude components for (intact) faces with different durations and high confidence (left hemisphere):
contra_high_N170_spec <- contrast(emmeans(lmm.opt.N170.high.E205, ~ Category + DuraBlock | Hemisphere), #
interaction = "pairwise", adjust = "Bonferroni", infer = TRUE)
contra_high_N170_spec[1:3]
## Category_pairwise DuraBlock_pairwise Hemisphere estimate SE df lower.CL upper.CL t.ratio p.value
## face - house 17_1 - 17_2 Left -0.3448969 0.2869752 18209.03 -1.0319733 0.3421795 -1.202 0.6883
## face - house 17_1 - 200_2 Left 1.5555782 0.2363195 7112.00 0.9897001 2.1214563 6.583 <.0001
## face - house 17_2 - 200_2 Left 1.9004751 0.2846656 13731.61 1.2189079 2.5820424 6.676 <.0001
##
## Degrees-of-freedom method: satterthwaite
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 3 estimates
## P value adjustment: bonferroni method for 3 tests
summary(contra_high_N170_spec[1:3], delta = delta_null, level = .9)
summary(contra_high_N170_spec[1:3], delta = delta_null, side = ">")
summary(contra_high_N170_spec[1:3], delta = delta_null, side = "<")
Comparisons of face-specific N170 amplitude components for (intact) faces with different durations and high confidence (right hemisphere):
# right
contra_high_N170_spec[4:6]
## Category_pairwise DuraBlock_pairwise Hemisphere estimate SE df lower.CL upper.CL t.ratio p.value
## face - house 17_1 - 17_2 Right 0.0542572 0.2862422 17571.01 -0.6310664 0.7395809 0.190 1.0000
## face - house 17_1 - 200_2 Right 1.9700598 0.2297485 2326.23 1.4196483 2.5204714 8.575 <.0001
## face - house 17_2 - 200_2 Right 1.9158026 0.2809292 7615.60 1.2431151 2.5884901 6.820 <.0001
##
## Degrees-of-freedom method: satterthwaite
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 3 estimates
## P value adjustment: bonferroni method for 3 tests
summary(contra_high_N170_spec[4:6], delta = delta_null, level = .9)
summary(contra_high_N170_spec[4:6], delta = delta_null, side = ">")
summary(contra_high_N170_spec[4:6], delta = delta_null, side = "<")
Comparisons of face-specific N170 amplitude component lateralization for (intact) faces with different durations and high confidence:
contra_high_N170_spec_hemi <- contrast(emmeans(lmm.opt.N170.high.E205, ~ Category + DuraBlock + Hemisphere), #
interaction = "pairwise", adjust = "Bonferroni", infer = TRUE)
contra_high_N170_spec_hemi
## Category_pairwise DuraBlock_pairwise Hemisphere_pairwise estimate SE df lower.CL upper.CL t.ratio p.value
## face - house 17_1 - 17_2 Left - Right -0.3991541 0.4044924 17388.02 -1.3675945 0.5692862 -0.987 0.9713
## face - house 17_1 - 200_2 Left - Right -0.4144816 0.3231755 2014.65 -1.1888040 0.3598408 -1.283 0.5994
## face - house 17_2 - 200_2 Left - Right -0.0153275 0.3961737 6959.66 -0.9639887 0.9333338 -0.039 1.0000
##
## Degrees-of-freedom method: satterthwaite
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 3 estimates
## P value adjustment: bonferroni method for 3 tests
summary(contra_high_N170_spec_hemi, delta = delta_null, level = .9)
summary(contra_high_N170_spec_hemi, delta = delta_null, side = ">")
summary(contra_high_N170_spec_hemi, delta = delta_null, side = "<")
Other results
Comparisons of P1 amplitudes for different durations and high confidence (without multiple comparison corrections):
summary(contra_high_N170[1:12], adjust = "none")
summary(contra_high_N170[1:12], delta = delta_null, adjust = "none", level = .9)
contra_high_N170_spec.hemi <- contrast(emmeans(lmm.opt.N170.high.E205, ~ Category + DuraBlock + Hemisphere), #
interaction = "pairwise", adjust = "Bonferroni", infer = TRUE)
contra_high_N170_spec.hemi
## Category_pairwise DuraBlock_pairwise Hemisphere_pairwise estimate SE df lower.CL upper.CL t.ratio p.value
## face - house 17_1 - 17_2 Left - Right -0.3991541 0.4044924 17388.02 -1.3675945 0.5692862 -0.987 0.9713
## face - house 17_1 - 200_2 Left - Right -0.4144816 0.3231755 2014.65 -1.1888040 0.3598408 -1.283 0.5994
## face - house 17_2 - 200_2 Left - Right -0.0153275 0.3961737 6959.66 -0.9639887 0.9333338 -0.039 1.0000
##
## Degrees-of-freedom method: satterthwaite
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 3 estimates
## P value adjustment: bonferroni method for 3 tests
summary(contra_high_N170_spec.hemi, delta = delta_null, level = .9)
Appendix
Model selection procedure
For the analysis for every dependent variable, an optimal model was obtained by the following general steps:
- Maximum Model: a maximum model with all fixed and random factors was built.
- ZCP Model: a zero-correlated-parameter (zcp) model based on the maximum model is built. The only difference between zcp and maximum models is that the correlations between random effects are forced to be zero in the zcp model.
- Reduced Model: the random effects which are not supported by the data will be removed from the zcp model. The function
step(fit, fixed.reduce = FALSE) from library(lmerTest) is used for this step. The algorithm could be found here.
- Extended Model: extending the reduced model with correlation parameters between the remaining random effects.
More details about this whole process could be found here: Bates, D., Kliegl, R., Vasishth, S., & Baayen, H. (2015). Parsimonious Mixed Models. Retrieved from http://arxiv.org/abs/1506.04967.
Stimuli specific effects
# 28 participants in total
resp.stim <- {
clean_beha_205 %>%
# filter(Duration == "17", Resp == "1") %>%
mutate(isKey1 = as.numeric(Resp == "1"),
Type_Category = paste(Type, Category, sep = "-")) %>%
group_by(Stimuli, Duration, Type_Category) %>%
summarize(avgCount = mean(isKey1))
}
ggplot(resp.stim, aes(x = Stimuli, y = avgCount, fill = Type_Category)) +
geom_col() +
facet_grid(Duration ~ .) +
labs(title = "Averaged percentage of being responsed as 'Key 1' for each stimuli",
y = "Percentage") +
theme(axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
text = element_text(size = 10),
axis.title = element_text(size = 16),
axis.text = element_text(size = 15), # the size of the texts in plot
# axis.text.x = element_text(angle = 45, vjust = 0.5),
legend.title=element_text(size=16),
legend.text=element_text(size=16),
# legend.position = "bottom",
legend.key.width = unit(1.2, "cm"),
plot.title = element_text(lineheight=.8, face="bold", size = 17),
strip.text = element_text(face="bold", size=15, lineheight=5.0),
strip.background = element_rect(fill="white", colour="white", size=1),
strip.placement = "outside",
legend.position = "bottom")
resp.stim.1 <- {
clean_beha_205 %>%
mutate(isKey1 = as.numeric(Resp == "1")) %>%
filter(Type == "intact", Category == "face") %>%
group_by(Stimuli, SubjCode, Duration) %>%
summarize(avgCount = mean(isKey1)) %>%
ungroup()
}
resp.stim.1.17 <- {
resp.stim.1 %>%
filter(Duration == "17") %>%
select(-Duration) %>%
spread(Stimuli, avgCount)
}
resp.stim.1.17 <- data.matrix(select(resp.stim.1.17, -SubjCode))
resp.stim.1.200 <- {
resp.stim.1 %>%
filter(Duration == "200") %>%
select(-Duration) %>%
spread(Stimuli, avgCount)
}
resp.stim.1.200 <- data.matrix(select(resp.stim.1.200, -SubjCode))
library(plotly)
heatmap_17 <- plot_ly(z = data.matrix(resp.stim.1.17), type = "heatmap", name = "Figrename")
heatmap_200 <- plot_ly(z = data.matrix(resp.stim.1.200), type = "heatmap", name = "Figrename")
Gray plots for publication
Equivalence tests
gg_equi_bw
Behavioural responses
plot_two_response
P1 – all
p1_all_LinePlot_bw <- {
data.frame(emm.P1.E205) %>%
mutate(Duration = fct_recode(Duration, `33` = "17", `216` = "200")) %>%
ggplot(aes(y = emmean, x = Duration, group = Category)) + # set the data, varialbes for x and y axises, and the variable for dividing data
geom_point(aes(fill = Duration), size=4, shape=21, stroke=0) + # , alpha = .7
scale_fill_manual(values = c("gray70", "black"), guide = "none") +
geom_line(aes(linetype = Category), size = .85) + # position = "dodge", alpha = .7
scale_linetype_manual(values=c("solid", "dashed")) +
geom_errorbar(mapping = aes(ymin = lower.CL, ymax = upper.CL), linetype = 1, # set the error bar
show.legend = FALSE, width = 0.15, alpha = .5) + # , position = position_dodge(width=0.9)
facet_grid(Type ~ Hemisphere, # switch = "x",
labeller = labeller(Hemisphere = hemiLabel)) +
scale_y_reverse(limits = ylimit_p1) + # set the limit for y axis
# scale_y_continuous(expand= c(0, 0)) + # remove the space between columns and x axis
labs(x = "Duration (ms)", y = expression(paste("P1 Amplitudes (", mu, "V)")), color = "Category", linetype = "Category") + # set the names for main, x and y axises title = "",
# scale_x_discrete(labels = c("left", "right")) +
# geom_text(label = c("", "", "", "", "***", "", "", ""), color = "red", size = 6, nudge_y = 0.5, nudge_x = 0.5) + # add starts to the significant columns
theme_bw() +
erp_theme +
theme(
legend.position = c(.5, .5),
legend.title=element_text(size=10),
legend.text=element_text(size=10),
) +
guides(linetype = guide_legend(title.position = "left", nrow = 1),
color = guide_legend(title.position = "left", nrow = 1)) +
NULL
}
p1_all_LinePlot_bw
N170 – all
n170_all_LinePlot_bw <- {
data.frame(emm.N170.E205) %>%
mutate(Duration = fct_recode(Duration, `33` = "17", `216` = "200")) %>%
ggplot(data = , aes(y = emmean, x = Duration, group = Category)) + # set the data, varialbes for x and y axises, and the variable for dividing data
geom_point(aes(fill = Duration), size=4, shape=21, stroke=0) + # , alpha = .7
scale_fill_manual(values = c("gray70", "black"), guide = "none") +
geom_line(aes(linetype = Category), size = .85) + # position = "dodge", alpha = .7
scale_linetype_manual(values=c("solid", "dashed")) +
geom_errorbar(mapping = aes(ymin = lower.CL, ymax = upper.CL), linetype = 1, # set the error bar
show.legend = FALSE, width = 0.15, alpha = .5) + # , position = position_dodge(width=0.9)
facet_grid(Type ~ Hemisphere, # switch = "x",
labeller = labeller(Hemisphere = hemiLabel)) +
scale_y_reverse(limits = ylimit_n170, breaks = seq(1, -5, -1)) + # set the limit for y axis
# scale_y_continuous(expand= c(0, 0)) + # remove the space between columns and x axis
labs(x = "Duration (ms)", y = expression(paste("N170 Amplitudes (", mu, "V)")), color = "Category", linetype = "Category") + # set the names for main, x and y axises title = "",
# scale_x_discrete(labels = c("left", "right")) +
# geom_text(label = c("", "", "", "", "***", "", "", ""), color = "red", size = 6, nudge_y = 0.5, nudge_x = 0.5) + # add starts to the significant columns
theme_bw() +
erp_theme +
theme(
legend.position = c(.5, .5),
legend.title=element_text(size=10),
legend.text=element_text(size=10),
) +
guides(linetype = guide_legend(title.position = "left", nrow = 1),
color = guide_legend(title.position = "left", nrow = 1))
}
n170_all_LinePlot_bw
P1 – con
p1_con_LinePlot_bw <- {
data.frame(emm.P1.SC.E205) %>%
mutate(DuraConf = fct_recode(DuraConf,
`33_guess` = "17_guess", `33_low` = "17_low", `33_high` = "17_high", `216_high` = "200_high")) %>%
ggplot(aes(y = emmean, x = DuraConf, group = Category)) + # set the data, varialbes for x and y axises, and the variable for dividing data
geom_point(aes(fill = DuraConf), size=4, shape=21, stroke=0) +
scale_fill_manual(values = c("gray75", "gray60", "gray35", "black"), guide = "none") +
geom_line(aes(linetype = Category), size = .85) + # position = "dodge", alpha = .7
scale_linetype_manual(values=c("solid", "dashed")) +
geom_errorbar(mapping = aes(ymin = lower.CL, ymax = upper.CL), linetype = 1, # set the error bar
show.legend = FALSE, width = 0.15, alpha = .5) + # , position = position_dodge(width=0.9)
facet_grid(. ~ Hemisphere, # switch = "x",
labeller = labeller(Hemisphere = hemiLabel)) +
scale_y_reverse(limits = ylimit_p1) + # set the limit for y axis
scale_x_discrete(labels = c("33 ms\nguess", "33 ms\nlow", "33 ms\nhigh", "216 ms\nhigh")) +
# scale_y_continuous(expand= c(0, 0)) + # remove the space between columns and x axis
labs(x = "Duration and Subjective Confidence", y = expression(paste("P1 Amplitudes (", mu, "V)")), fill = "Category") + # set the names for main, x and y axises title = "",
theme_bw() +
erp_theme +
theme(
axis.title.x = element_text(margin = margin(t = 10, r = 0, b = 0, l = 0)),
legend.margin = margin(t = -10),
legend.title=element_text(size=12),
legend.text=element_text(size=12),
)
}
p1_con_LinePlot_bw
N170 – con
n170_con_LinePlot_bw = {
data.frame(emm.N170.SC.E205) %>%
mutate(DuraConf = fct_recode(DuraConf,
`33_guess` = "17_guess", `33_low` = "17_low", `33_high` = "17_high", `216_high` = "200_high")) %>%
ggplot(aes(y = emmean, x = DuraConf, group = Category)) + # set the data, varialbes for x and y axises, and the variable for dividing data
geom_point(aes(fill = DuraConf), size=4, shape=21, stroke=0) +
scale_fill_manual(values = c("gray75", "gray60", "gray35", "black"), guide = "none") +
geom_line(aes(linetype = Category), size = .85) + # position = "dodge", alpha = .7
scale_linetype_manual(values=c("solid", "dashed")) +
geom_errorbar(mapping = aes(ymin = lower.CL, ymax = upper.CL), linetype = 1, # set the error bar
show.legend = FALSE, width = 0.15, alpha = .5) + # , position = position_dodge(width=0.9)
facet_grid(. ~ Hemisphere, # switch = "x",
labeller = labeller(Hemisphere = hemiLabel)) +
scale_y_reverse(limits = ylimit_n170, breaks = seq(1, -5, -1)) + # set the limit for y axis
scale_x_discrete(labels = c("33 ms\nguess", "33 ms\nlow", "33 ms\nhigh", "216 ms\nhigh")) +
# scale_y_continuous(expand= c(0, 0)) + # remove the space between columns and x axis
labs(x = "Duration and Subjective Confidence", y = expression(paste("N170 Amplitudes (", mu, "V)")), fill = "Category") + # set the names for main, x and y axises title = "",
theme_bw() +
erp_theme +
theme(
axis.title.x = element_text(margin = margin(t = 10, r = 0, b = 0, l = 0)),
legend.margin = margin(t = -10),
legend.title=element_text(size=12),
legend.text=element_text(size=12),
)
}
n170_con_LinePlot_bw
P1 – high
p1_high_LinePlot_bw = {
data.frame(emm.P1.high.E205) %>%
mutate(DuraBlock = fct_recode(DuraBlock,
`33_1` = "17_1", `33_2` = "17_2", `216_2` = "200_2")) %>%
ggplot(data = , aes(y = emmean, x = DuraBlock, group = Category)) + # set the data, variables for x and y axises, and the variable for dividing data
geom_point(aes(fill = DuraBlock), size=4, shape=21, stroke=0) +
scale_fill_manual(values = c("gray60", "gray35", "black"), guide = "none") +
geom_line(aes(linetype = Category), size = .85) + # position = "dodge", alpha = .7
scale_linetype_manual(values=c("solid", "dashed")) +
geom_errorbar(mapping = aes(ymin = lower.CL, ymax = upper.CL), linetype = 1, # set the error bar
show.legend = FALSE, width = 0.15, alpha = .5) + # , position = position_dodge(width=0.9)
facet_grid(. ~ Hemisphere, # switch = "x",
labeller = labeller(Hemisphere = hemiLabel)) +
scale_y_reverse(limits = ylimit_p1) + # set the limit for y axis
scale_x_discrete(labels = c("33 ms\nhalf 1", "33 ms\nhalf 2", "216 ms\nhalf 2")) +
# scale_y_continuous(expand= c(0, 0)) + # remove the space between columns and x axis
labs(x = "Duration and Halves", y = expression(paste("P1 Amplitudes (", mu, "V)")), fill = "Category") + # set the names for main, x and y axises title = "",
theme_bw() +
erp_theme +
theme(
axis.title.x = element_text(margin = margin(t = 10, r = 0, b = 0, l = 0)),
legend.margin = margin(t = -10),
legend.title=element_text(size=12),
legend.text=element_text(size=12),
)
}
p1_high_LinePlot_bw
N170 – high
n170_high_LinePlot_bw = {
data.frame(emm.N170.high.E205) %>%
mutate(DuraBlock = fct_recode(DuraBlock,
`33_1` = "17_1", `33_2` = "17_2", `216_2` = "200_2")) %>%
ggplot(aes(y = emmean, x = DuraBlock, group = Category)) + # set the data, varialbes for x and y axises, and the variable for dividing data
geom_point(aes(fill = DuraBlock), size=4, shape=21, stroke=0) +
scale_fill_manual(values = c("gray60", "gray35", "black"), guide = "none") +
geom_line(aes(linetype = Category), size = .85) + # position = "dodge", alpha = .7
scale_linetype_manual(values=c("solid", "dashed")) +
geom_errorbar(mapping = aes(ymin = lower.CL, ymax = upper.CL), linetype = 1, # set the error bar
show.legend = FALSE, width = 0.15, alpha = .5) + # , position = position_dodge(width=0.9)
facet_grid(. ~ Hemisphere, # switch = "x",
labeller = labeller(Hemisphere = hemiLabel)) +
scale_y_reverse(limits = ylimit_n170, breaks = seq(1, -5, -1)) + # set the limit for y axis
scale_x_discrete(labels = c("33 ms\nhalf 1", "33 ms\nhalf 2", "216 ms\nhalf 2")) +
# scale_y_continuous(expand= c(0, 0)) + # remove the space between columns and x axis
labs(x = "Duration and Halves", y = expression(paste("N170 Amplitudes (", mu, "V)")), fill = "Category") + # set the names for main, x and y axises title = "",
theme_bw() +
erp_theme +
theme(
axis.title.x = element_text(margin = margin(t = 10, r = 0, b = 0, l = 0)),
legend.margin = margin(t = -10),
legend.title=element_text(size=12),
legend.text=element_text(size=12),
)
}
n170_high_LinePlot_bw
Save plots
# filetype <- "pdf"
# ggsave(paste0("equivalence_test_bw.", filetype), gg_equi_bw, width = 8, height = 4.7)
# ggsave(paste0('plot_two_response_bw.', filetype), plot_two_response, width = 10, height = 5)
# ggsave(paste0('P1_all_pub_bw.', filetype), p1_all_LinePlot_bw, width = 6, height = 5)
# ggsave(paste0('N170_all_pub_bw.', filetype), n170_all_LinePlot_bw, width = 6, height = 5)
# ggsave(paste0('P1_con_pub_bw.', filetype), p1_con_LinePlot_bw, width = 7, height = 5)
# ggsave(paste0('N170_con_pub_bw.', filetype), n170_con_LinePlot_bw, width = 7, height = 5)
# ggsave(paste0('p1_high_pub_bw.', filetype), p1_high_LinePlot_bw, width = 7, height = 4.5)
# ggsave(paste0('n170_high_pub_bw.', filetype), n170_high_LinePlot_bw, width = 7, height = 4.5)
Stimulus duration
The duration in this section was the differences between the onset of masking and the onset of experimental stimuli, which were the time points of triggers sent by E-prime to EGI amplifier.
stim_dura <- read_csv(file.path("data", "E205_stimulus_duration.csv"))
stim_dura %>%
group_by(OnEvent, ISI) %>%
summarize(count = n())
As shown in the above table, when the stimuli were intended to presented for 17ms, the actual interval between the onset of stimulus and masking were 32 (only a few trials), 33, or 34ms. When the stimuli were intended to presented for 200ms, the actual interval between the onset of stimulus and masking were 216 or 217 ms.
Timing test
The
timing <- read_csv(file.path("data", "timing_test_duration.csv"))
timing %>%
group_by(OnEvents, timing) %>%
summarize(count = n())
timing %>%
group_by(OnEvents, dura) %>%
summarize(count = n())
Reply to Reviewer
# ("acc") descriptive statistics of accuracy for plotting
acc.dist.205 <- clean_beha_205 %>%
mutate(Resp_acc = if_else(Resp %in% c(1,2), "1 or 2",
if_else(Resp %in% c(4,5), "4 or 5",
if_else(Resp == 3, "3", "NA")))) %>%
group_by(SubjCode, Type, Category, Duration, Resp_acc) %>%
summarize(Count = n(), .groups = "drop") %>%
right_join(., tidyr::expand(., SubjCode, Type, Category, Duration, Resp_acc)) %>%
mutate(Count = if_else(is.na(Count), 0L, Count)) %>%
right_join(tmp.count, by = c("SubjCode", "Type", "Category", "Duration")) %>%
mutate(RespRate = Count / Count_sum) %>%
# ungroup()
group_by(Type, Category, Duration, Resp_acc) %>%
summarize(N = n()
, Mean = sum(RespRate)/N
, SD = sd(RespRate)
, SE = SD/sqrt(N)
, SE.lo = Mean - SE, SE.hi = Mean + SE
, CI_lo = Mean - SE * 1.96
, CI_hi = Mean + SE * 1.96
, Median = median(RespRate)
, Lower = Mean - SD, Upper = Mean + SD # for rainclound plot
)
acc.dist.205
dura.labels <- c("17ms", "200ms")
names(dura.labels) <- c(17, 200)
dist.acc.ColuPlot.205 <- {
ggplot(data = acc.dist.205, aes(y = Mean, x = Resp_acc, fill = Category)) + # set the data, varialbes for x and y axises, and the variable for dividing data
geom_col(position = "dodge", color = "black", alpha = .7) + # position of columns and countour of columns
facet_grid(Type ~ Duration,
labeller = labeller(Duration = dura.labels)) +
geom_errorbar(mapping = aes(ymin = CI_lo, ymax = CI_hi), linetype = 1, # set the error bar
show.legend = FALSE, width = 0.25, alpha = .5,
position = position_dodge(width=0.9)) +
geom_hline(yintercept = c(1), linetype = 5, alpha = 0.5) + # add the line for 0.5 and 1 (y)
coord_cartesian(ylim = c(0,1.1)) + # set the limit for y axis
scale_y_continuous(expand= c(0, 0)) + # remove the space between columns and x axis
labs(title = "", x = "Responses", y = "Percentage", fill = "Category") + # set the names for main, x and y axises Subjective Responses for E205
scale_fill_grey() + # set the color for columns
# geom_text(label = c("***", "***", "", ""), color = c("red", "red", "red", "red"), size = 6, nudge_y = 0.05, nudge_x = rep(c(-0.225, 0.225), 2)) + # add starts to the significant columns
erp_theme
}
# ggsave("acc_column_plot.png", dist.acc.ColuPlot.205, width = 8, height = 7)
dist.acc.ColuPlot.205
