library(tidyverse)
## Warning: package 'tidyr' was built under R version 4.3.3
## Warning: package 'readr' was built under R version 4.3.3
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(bruceR)
##
## bruceR (v2023.9)
## Broadly Useful Convenient and Efficient R functions
##
## Packages also loaded:
## ✔ data.table ✔ emmeans
## ✔ dplyr ✔ lmerTest
## ✔ tidyr ✔ effectsize
## ✔ stringr ✔ performance
## ✔ ggplot2 ✔ interactions
##
## Main functions of `bruceR`:
## cc() Describe() TTEST()
## add() Freq() MANOVA()
## .mean() Corr() EMMEANS()
## set.wd() Alpha() PROCESS()
## import() EFA() model_summary()
## print_table() CFA() lavaan_summary()
##
## For full functionality, please install all dependencies:
## install.packages("bruceR", dep=TRUE)
##
## Online documentation:
## https://psychbruce.github.io/bruceR
##
## To use this package in publications, please cite:
## Bao, H.-W.-S. (2023). bruceR: Broadly useful convenient and efficient R functions (Version 2023.9) [Computer software]. https://CRAN.R-project.org/package=bruceR
##
##
## NEWS: A new version of bruceR (2024.6) is available (2024-06-13)!
##
## ***** Please update *****
## install.packages("bruceR", dep=TRUE)
##
##
## These packages are dependencies of `bruceR` but not installed:
## - pacman, ggtext, see, vars, phia, GGally, GPArotation
##
## ***** Install all dependencies *****
## install.packages("bruceR", dep=TRUE)
##
##
## Attaching package: 'bruceR'
##
## The following object is masked _by_ 'package:data.table':
##
## %notin%
#清空环境变量
rm(list = ls())
set.wd()
## ✔ Set working directory to "C:/Users/psyuser/Desktop/CP conference"
source('summarySE.R')
library(afex)
## Warning: package 'afex' was built under R version 4.3.3
## ************
## Welcome to afex. For support visit: http://afex.singmann.science/
## - Functions for ANOVAs: aov_car(), aov_ez(), and aov_4()
## - Methods for calculating p-values with mixed(): 'S', 'KR', 'LRT', and 'PB'
## - 'afex_aov' and 'mixed' objects can be passed to emmeans() for follow-up tests
## - Get and set global package options with: afex_options()
## - Set sum-to-zero contrasts globally: set_sum_contrasts()
## - For example analyses see: browseVignettes("afex")
## ************
##
## Attaching package: 'afex'
##
## The following object is masked from 'package:lme4':
##
## lmer
library(emmeans) # for post hoc test
library(ggthemes)
## Warning: package 'ggthemes' was built under R version 4.3.3
library(cowplot)
## Warning: package 'cowplot' was built under R version 4.3.3
##
## Attaching package: 'cowplot'
##
## The following object is masked from 'package:ggthemes':
##
## theme_map
##
## The following object is masked from 'package:lubridate':
##
## stamp
library(ggpubr)
## Registered S3 methods overwritten by 'broom':
## method from
## tidy.glht jtools
## tidy.summary.glht jtools
##
## Attaching package: 'ggpubr'
##
## The following object is masked from 'package:cowplot':
##
## get_legend
library(rstatix)
##
## Attaching package: 'rstatix'
##
## The following objects are masked from 'package:effectsize':
##
## cohens_d, eta_squared
##
## The following object is masked from 'package:stats':
##
## filter
library(effsize)
library(data.table)
my_data <- read_csv('edfData_DotFace_Allsubject_in_trial_fix.csv') %>% filter(subnum != 'HC029_')
## Rows: 136014 Columns: 29
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (6): subnum, sub_type, expression, L_stim, R_stim, emo_stim
## dbl (23): index, startT, endT, duration, avgX, avgY, avgPupil, AOI, AOI_gaze...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
my_data <- my_data %>% mutate(trial_type=case_when(
trial_type == "1" ~ "guided_emo",
trial_type == "2" ~ "guided_neu",
trial_type == "3" ~ "centered",
TRUE ~ NA_character_)) %>%
mutate(dwell_emo=case_when(
dwell_emo == "1" ~ "guided_emo",
dwell_emo == "2" ~ "guided_neu",
dwell_emo == "31" ~ "centered_emo",
dwell_emo == "32" ~ "centered_neu")) %>%
filter(duration >= 60) %>%
group_by(subnum,Ntrial) %>%
mutate(continue_gaze=rleid(AOI_emo)) %>%
ungroup()
# 看看有效引导的百分比
my_data_1ratio <- my_data %>% group_by(subnum, sub_type, trial_type) %>%
select(c(subnum, sub_type, trial_type, guided_trial)) %>%
replace_na(list(guided_trial = 0)) %>%
summarise(across(everything(), ~mean(. == 1))) %>% ungroup()
## `summarise()` has grouped output by 'subnum', 'sub_type'. You can override
## using the `.groups` argument.
my_data <- my_data %>% filter(subnum != 'MD034_') # filtered because the guided trial rate is lower than 70%
# 首注视偏好
# 每个被试每种表情的注视偏好比例
data4 <- my_data %>%
filter(trial_type == 'centered', AOI_emo==1|AOI_emo==2) %>%
group_by(subnum, Ntrial, sub_type, expression, AOI_emo) %>%
mutate(order=1:n()) %>%
ungroup() %>%
filter(order==1) %>%
group_by(subnum, Ntrial, sub_type, expression) %>%
mutate(order=1:n(),AOI_emo=factor(AOI_emo,levels=c(1,2),
labels=c('emo','neu'))) %>%
ungroup() %>%
filter(order==1) %>%
group_by(subnum, sub_type, expression, AOI_emo) %>%
dplyr::summarize(num=n()) %>%
ungroup() %>%
spread(key = AOI_emo, value = num) %>%
mutate(tot=emo+neu,rate=emo/tot)
## `summarise()` has grouped output by 'subnum', 'sub_type', 'expression'. You can
## override using the `.groups` argument.
# Test if the first gaze rate within each bs and ws factor is significantly higher than 0.5
t_test_results <- data4 %>%
group_by(sub_type, expression) %>%
summarise(t_test = list(t.test(rate, mu = 0.5, alternative = "greater"))) %>%
mutate(p_value = map_dbl(t_test, "p.value"),
statistic = map_dbl(t_test, "statistic"),
conf_int = map(t_test, "conf.int"))
## `summarise()` has grouped output by 'sub_type'. You can override using the
## `.groups` argument.
# Apply Bonferroni correction
num_tests <- nrow(t_test_results)
t_test_results <- t_test_results %>%
mutate(p_value_bonferroni = p.adjust(p_value, method = "bonferroni", n = num_tests))
# Print the t-test results
print(t_test_results)
## # A tibble: 6 × 7
## # Groups: sub_type [2]
## sub_type expression t_test p_value statistic conf_int p_value_bonferroni
## <chr> <chr> <list> <dbl> <dbl> <list> <dbl>
## 1 HC F <htest> 3.07e- 8 6.35 <dbl [2]> 1.84e- 7
## 2 HC H <htest> 3.33e- 1 0.435 <dbl [2]> 1 e+ 0
## 3 HC S <htest> 4.83e- 3 2.69 <dbl [2]> 2.90e- 2
## 4 MDD F <htest> 8.42e-11 7.41 <dbl [2]> 5.05e-10
## 5 MDD H <htest> 1.36e- 3 3.10 <dbl [2]> 8.14e- 3
## 6 MDD S <htest> 4.31e- 8 5.94 <dbl [2]> 2.58e- 7
# Perform ANOVA for first fixation preference
anova_first_prefer <- aov_ez(
data = data4,
id = 'subnum',
dv = 'rate',
between = 'sub_type',
within = 'expression',
observed = 'sub_type'
)
## Converting to factor: sub_type
## Contrasts set to contr.sum for the following variables: sub_type
# Summarize the ANOVA results
summary(anova_first_prefer)
## Warning in summary.Anova.mlm(object$Anova, multivariate = FALSE): HF eps > 1
## treated as 1
##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 116.326 1 2.0477 124 7044.1570 < 2.2e-16 ***
## sub_type 0.029 1 2.0477 124 1.7853 0.1840
## expression 0.365 2 3.5539 248 12.7495 5.364e-06 ***
## sub_type:expression 0.032 2 3.5539 248 1.1222 0.3272
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## expression 0.99883 0.93029
## sub_type:expression 0.99883 0.93029
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## expression 0.99883 5.421e-06 ***
## sub_type:expression 0.99883 0.3272
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## expression 1.015172 5.364301e-06
## sub_type:expression 1.015172 3.271978e-01
# Perform post hoc test on 'expression' variable
emmeans_result <- emmeans(anova_first_prefer, pairwise ~ expression, adjust = "bonferroni")
# Print the post hoc test results
print(emmeans_result)
## $emmeans
## expression emmean SE df lower.CL upper.CL
## F 0.604 0.0108 124 0.583 0.626
## H 0.527 0.0114 124 0.504 0.549
## S 0.564 0.0112 124 0.542 0.587
##
## Results are averaged over the levels of: sub_type
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## F - H 0.0776 0.0151 124 5.138 <.0001
## F - S 0.0397 0.0155 124 2.566 0.0344
## H - S -0.0379 0.0155 124 -2.441 0.0482
##
## Results are averaged over the levels of: sub_type
## P value adjustment: bonferroni method for 3 tests
# 引导中间 潜伏期
outlier.IQR <- function(x, multiple = 1.5) {
q <- quantile(x, na.rm = TRUE) #四分位间距3倍间距以外的认为是离群值
IQR <- q[4] - q[2]
x1 <- q[2] - multiple * IQR
x2 <- q[4] + multiple * IQR
return(c(x1, x2))
}
data3 <- my_data %>%
filter(trial_type == 'centered', AOI_emo==1|AOI_emo==2) %>%
group_by(subnum, Ntrial, sub_type, expression, AOI_emo) %>%
mutate(order=1:n()) %>%
ungroup() %>%
filter(order==1) %>%
mutate(latency=startT-Tstart,AOI_emo=factor(AOI_emo)) %>%
# filter(latency<=2500, latency>=120) %>% # need papers to support
group_by(subnum, Ntrial, sub_type, trial_type, expression) %>%
summarise(
emo_latency = latency[AOI_emo == 1],
neu_latency = latency[AOI_emo == 2],
emo_neu_latency = emo_latency - neu_latency
) %>%
ungroup() %>%
# kick outlier within subject
group_by(subnum) %>%
mutate(
IQR_L = outlier.IQR(emo_neu_latency)[1],
IQR_H = outlier.IQR(emo_neu_latency)[2],
valid1 = if_else(!is.na(IQR_L) & !is.na(IQR_H) & emo_neu_latency >= IQR_L & emo_neu_latency <= IQR_H, 1, 0)
) %>%
filter(valid1 == 1) %>%
ungroup()
## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## `summarise()` has grouped output by 'subnum', 'Ntrial', 'sub_type',
## 'trial_type', 'expression'. You can override using the `.groups` argument.
sumrepdat3E_sub <- summarySE(data3, measurevar = "emo_neu_latency",
groupvars=c('subnum',"sub_type","trial_type","expression"))
## Warning in qt(conf.interval/2 + 0.5, datac$N - 1): NaNs produced
# Test if the emo - neu first fixation latency within each bs and ws factor is significantly higher than 0.5
t_test_results <- data3 %>%
group_by(sub_type, expression) %>%
summarise(t_test = list(t.test(emo_neu_latency, mu = 0, alternative = "two.sided"))) %>%
mutate(p_value = map_dbl(t_test, "p.value"),
statistic = map_dbl(t_test, "statistic"),
conf_int = map(t_test, "conf.int"))
## `summarise()` has grouped output by 'sub_type'. You can override using the
## `.groups` argument.
# Apply Bonferroni correction
num_tests <- nrow(t_test_results)
t_test_results <- t_test_results %>%
mutate(p_value_bonferroni = p.adjust(p_value, method = "bonferroni", n = num_tests))
# Print the t-test results
print(t_test_results)
## # A tibble: 6 × 7
## # Groups: sub_type [2]
## sub_type expression t_test p_value statistic conf_int p_value_bonferroni
## <chr> <chr> <list> <dbl> <dbl> <list> <dbl>
## 1 HC F <htest> 6.69e-13 -7.32 <dbl [2]> 4.01e-12
## 2 HC H <htest> 9.49e- 2 -1.67 <dbl [2]> 5.70e- 1
## 3 HC S <htest> 1.75e- 4 -3.77 <dbl [2]> 1.05e- 3
## 4 MDD F <htest> 1.29e-15 -8.15 <dbl [2]> 7.76e-15
## 5 MDD H <htest> 3.20e- 4 -3.61 <dbl [2]> 1.92e- 3
## 6 MDD S <htest> 7.95e- 5 -3.96 <dbl [2]> 4.77e- 4
# Perform ANOVA for first fixation latency
anova_fix_latency <- aov_ez(
data = sumrepdat3E_sub,
id = 'subnum',
dv = 'emo_neu_latency_mean',
between = 'sub_type',
within = 'expression',
observed = 'sub_type'
)
## Converting to factor: sub_type
## Warning: Missing values for 2 ID(s), which were removed before analysis:
## MD029_, MD039_
## Below the first few rows (in wide format) of the removed cases with missing data.
## subnum sub_type F H S
## # 60 MD029_ MDD 220 NA NA
## # 64 MD039_ MDD NA -370.6667 -492
## Contrasts set to contr.sum for the following variables: sub_type
# Summarize the ANOVA results
summary(anova_fix_latency)
##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 5332435 1 7156544 122 90.9038 < 2.2e-16 ***
## sub_type 5161 1 7156544 122 0.0880 0.7672594
## expression 938153 2 12865476 244 8.8963 0.0001865 ***
## sub_type:expression 63988 2 12865476 244 0.6068 0.5459228
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## expression 0.96949 0.15346
## sub_type:expression 0.96949 0.15346
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## expression 0.9704 0.0002214 ***
## sub_type:expression 0.9704 0.5411534
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## expression 0.9858375 0.0002024835
## sub_type:expression 0.9858375 0.5436580898
# Perform post hoc test on 'sub_type' variable
emmeans_latency_result <- emmeans(anova_fix_latency, pairwise ~ expression, adjust = "bonferroni")
# Print the post hoc test results
print(emmeans_latency_result)
## $emmeans
## expression emmean SE df lower.CL upper.CL
## F -192.0 22.1 122 -236 -148.4
## H -72.7 22.2 122 -117 -28.7
## S -100.2 19.6 122 -139 -61.4
##
## Results are averaged over the levels of: sub_type
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## F - H -119.4 31.2 122 -3.824 0.0006
## F - S -91.8 30.5 122 -3.005 0.0097
## H - S 27.5 27.0 122 1.022 0.9263
##
## Results are averaged over the levels of: sub_type
## P value adjustment: bonferroni method for 3 tests
# 引导条件下首滞留
outlier.IQR <- function(x, multiple = 1.5) {
q <- quantile(x, na.rm = TRUE) #四分位间距3倍间距以外的认为是离群值
IQR <- q[4] - q[2]
x1 <- q[2] - multiple * IQR
x2 <- q[4] + multiple * IQR
return(c(x1, x2))
}
# Function to check continuity of index numbers within each group
check_continuity <- function(df) {
smallest_index <- min(df$index)
df <- df %>% arrange(index)
for (i in seq_along(df$index)) {
if (i == 1) next
if (df$index[i] != df$index[i - 1] + 1) {
return(df[1:(i - 1), ])
}
}
return(df)
}
# Function to process data for a given trial type and AOI_emo
process_data <- function(data, Trial_type, aoi_emo) {
data <- data %>%
filter(trial_type == Trial_type) %>%
filter(AOI_emo == aoi_emo) %>%
group_by(subnum, Ntrial) %>%
group_modify(~ check_continuity(.x)) %>%
mutate(dwell_time = sum(duration, na.rm = TRUE)) %>%
ungroup()
return(data)
}
# Process data for data2_emo
data2_emo <- process_data(my_data, 'guided_emo', 1)
# Create new variable data2_emo_dwell
data2_emo_dwell <- data2_emo %>%
select(subnum, sub_type, trial_type, Ntrial, expression, duration) %>%
group_by(subnum, Ntrial) %>%
mutate(dwell_emo = sum(duration, na.rm = TRUE)) %>%
ungroup() %>%
select(-duration) %>% # Drop the duration column
distinct(subnum, Ntrial, .keep_all = TRUE) # Drop duplicate rows based on subnum and Ntrial
# The first dwell time for guided_emo stimulus (emo dwell)
data2_emo_dwell_subed <- data2_emo_dwell %>%
# kick outlier within subject
group_by(subnum) %>%
mutate(
IQR_L = outlier.IQR(dwell_emo)[1],
IQR_H = outlier.IQR(dwell_emo)[2],
valid1 = if_else(!is.na(IQR_L) & !is.na(IQR_H) & dwell_emo >= IQR_L & dwell_emo <= IQR_H, 1, 0)
) %>%
filter(valid1 == 1) %>%
ungroup()
# Process data for data2_neu
data2_neu <- process_data(my_data, 'guided_neu', 2)
# Create new variable data2_neu_dwell
data2_neu_dwell <- data2_neu %>%
select(subnum, sub_type, trial_type, Ntrial, expression, duration) %>%
group_by(subnum, Ntrial) %>%
mutate(dwell_neu = sum(duration, na.rm = TRUE)) %>%
ungroup() %>%
select(-duration) %>% # Drop the duration column
distinct(subnum, Ntrial, .keep_all = TRUE) # Drop duplicate rows based on subnum and Ntrial
# The first dwell time for guided_neu stimulus (neu dwell)
data2_neu_dwell_subed <- data2_neu_dwell %>%
# kick outlier within subject
group_by(subnum) %>%
mutate(
IQR_L = outlier.IQR(dwell_neu)[1],
IQR_H = outlier.IQR(dwell_neu)[2],
valid1 = if_else(!is.na(IQR_L) & !is.na(IQR_H) & dwell_neu >= IQR_L & dwell_neu <= IQR_H, 1, 0)
) %>%
filter(valid1 == 1) %>%
ungroup()
# Summary statistics for emo_neu_dwell
sumrepdat2E_sub <- summarySE(data2_emo_dwell_subed, measurevar = "dwell_emo",
groupvars=c('subnum',"sub_type", "trial_type", "expression"))
sumrepdat2N_sub <- summarySE(data2_neu_dwell_subed, measurevar = "dwell_neu",
groupvars=c('subnum',"sub_type", "trial_type", "expression"))
# Perform ANOVA for first dwell on emo
anova_emo_dwell <- aov_ez(
data = sumrepdat2E_sub,
id = 'subnum',
dv = 'dwell_emo_mean',
between = 'sub_type',
within = 'expression',
observed = 'sub_type'
)
## Converting to factor: sub_type
## Contrasts set to contr.sum for the following variables: sub_type
# Summarize the ANOVA results
summary(anova_emo_dwell)
## Warning in summary.Anova.mlm(object$Anova, multivariate = FALSE): HF eps > 1
## treated as 1
##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 112743803 1 28280156 124 494.3477 < 2e-16 ***
## sub_type 315960 1 28280156 124 1.3854 0.24144
## expression 62413 2 3213528 248 2.4083 0.09207 .
## sub_type:expression 2336 2 3213528 248 0.0901 0.91382
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## expression 0.99308 0.65251
## sub_type:expression 0.99308 0.65251
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## expression 0.99313 0.09247 .
## sub_type:expression 0.99313 0.91270
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## expression 1.009241 0.09206763
## sub_type:expression 1.009241 0.91382407
# Perform post hoc test on 'expression' variable
emmeans_result <- emmeans(anova_emo_dwell, pairwise ~ expression, adjust = "bonferroni")
# Print the post hoc test results
print(emmeans_result)
## $emmeans
## expression emmean SE df lower.CL upper.CL
## F 567 26.5 124 515 620
## H 564 28.9 124 506 621
## S 538 23.6 124 491 585
##
## Results are averaged over the levels of: sub_type
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## F - H 3.81 15.1 124 0.253 1.0000
## F - S 29.47 14.0 124 2.100 0.1132
## H - S 25.67 14.7 124 1.746 0.2497
##
## Results are averaged over the levels of: sub_type
## P value adjustment: bonferroni method for 3 tests
# Perform ANOVA for first dwell on neu
anova_neu_dwell <- aov_ez(
data = sumrepdat2N_sub,
id = 'subnum',
dv = 'dwell_neu_mean',
between = 'sub_type',
within = 'expression',
observed = 'sub_type'
)
## Converting to factor: sub_type
## Contrasts set to contr.sum for the following variables: sub_type
# Summarize the ANOVA results
summary(anova_neu_dwell)
##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 103210986 1 23144143 124 552.9763 < 2e-16 ***
## sub_type 367794 1 23144143 124 1.9705 0.16289
## expression 36701 2 1460157 248 3.1167 0.04604 *
## sub_type:expression 14383 2 1460157 248 1.2214 0.29657
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## expression 0.69728 2.3418e-10
## sub_type:expression 0.69728 2.3418e-10
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## expression 0.76762 0.05983 .
## sub_type:expression 0.76762 0.28929
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## expression 0.7753503 0.0593184
## sub_type:expression 0.7753503 0.2896198
# Perform post hoc test on 'expression' variable
emmeans_result <- emmeans(anova_neu_dwell, pairwise ~ expression, adjust = "bonferroni")
# Print the post hoc test results
print(emmeans_result)
## $emmeans
## expression emmean SE df lower.CL upper.CL
## F 546 26.6 124 493 599
## H 529 22.7 124 484 574
## S 522 20.3 124 482 562
##
## Results are averaged over the levels of: sub_type
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## F - H 17.31 9.63 124 1.797 0.2242
## F - S 23.78 12.05 124 1.974 0.1520
## H - S 6.47 7.28 124 0.889 1.0000
##
## Results are averaged over the levels of: sub_type
## P value adjustment: bonferroni method for 3 tests
library(BayesFactor)
## Loading required package: coda
## Warning in .recacheSubclasses(def@className, def, env): undefined subclass
## "ndiMatrix" of class "replValueSp"; definition not updated
## ************
## Welcome to BayesFactor 0.9.12-4.7. If you have questions, please contact Richard Morey (richarddmorey@gmail.com).
##
## Type BFManual() to open the manual.
## ************
# Prepare the data
first_prefer_data <- data4
# Ensure subnum is a factor
first_prefer_data$subnum <- as.factor(first_prefer_data$subnum)
first_prefer_data$sub_type <- as.factor(first_prefer_data$sub_type)
first_prefer_data$expression <- as.factor(first_prefer_data$expression)
# Perform Bayesian ANOVA
bf_anova <- anovaBF(rate ~ sub_type * expression + subnum, data = first_prefer_data, whichRandom = "subnum")
## Warning: data coerced from tibble to data frame
# Print the results
print(bf_anova)
## Bayes factor analysis
## --------------
## [1] sub_type + subnum : 0.2831596 ±4.82%
## [2] expression + subnum : 2807.68 ±2.14%
## [3] sub_type + expression + subnum : 761.2163 ±1.91%
## [4] sub_type + expression + sub_type:expression + subnum : 118.1741 ±2.61%
##
## Against denominator:
## rate ~ subnum
## ---
## Bayes factor type: BFlinearModel, JZS
# Extract posterior samples from the full model (index 4)
# The full model includes sub_type, expression, their interaction, and subnum.
posterior_samples <- posterior(bf_anova, iterations = 1000, index = 4)
# # Inspect the column names in the posterior_samples object
# print(colnames(posterior_samples))
# Test specific hypothesis: HC in sub_type is different from MDD in sub_type
# Calculate the difference between HC and MDD in sub_type
diff_HC_MDD <- posterior_samples[, "sub_type-HC"] - posterior_samples[, "sub_type-MDD"]
# Calculate the Bayes factor for the hypothesis
bf_hypothesis <- mean(diff_HC_MDD > 0) / mean(diff_HC_MDD < 0)
# Print the hypothesis test results
print(bf_hypothesis) # 0.000, which suppport the hypothesis that HC in sub_type is different from MDD in sub_type
## [1] 0.1273957
#####
# Perform Bayesian ANOVA for emo_neu_latency
latency_data <- data3
# Ensure subnum is a factor
latency_data$subnum <- as.factor(latency_data$subnum)
latency_data$sub_type <- as.factor(latency_data$sub_type)
latency_data$expression <- as.factor(latency_data$expression)
# Perform Bayesian ANOVA
bf_anova_latency <- anovaBF(emo_neu_latency ~ sub_type * expression + subnum, data = latency_data, whichRandom = "subnum")
## Warning: data coerced from tibble to data frame
# Print the results
print(bf_anova_latency)
## Bayes factor analysis
## --------------
## [1] sub_type + subnum : 0.06182154 ±5.24%
## [2] expression + subnum : 715.5126 ±0.8%
## [3] sub_type + expression + subnum : 43.92082 ±1.74%
## [4] sub_type + expression + sub_type:expression + subnum : 0.3442246 ±1.67%
##
## Against denominator:
## emo_neu_latency ~ subnum
## ---
## Bayes factor type: BFlinearModel, JZS
# Extract posterior samples from the full model (index 4)
# The full model includes sub_type, expression, their interaction, and subnum.
posterior_samples_latency <- posterior(bf_anova_latency, iterations = 1000, index = 4)
# summary(posterior_samples_latency)
# sub_type-HC and sub_type-MDD:
# These terms represent the estimated effects (coefficients) for the levels HC and MDD of the sub_type factor, respectively.
# In a Bayesian ANOVA, each level of a factor (except the reference level) has its own coefficient that represents the effect of that level compared to the reference level.
# For example, if HC is the reference level, sub_type-MDD would represent the effect of being in the MDD group compared to the HC group.
# # Inspect the column names in the posterior_samples object
# print(colnames(posterior_samples_latency))
# Test specific hypothesis: HC in sub_type is different from MDD in sub_type
# Calculate the difference between HC and MDD in sub_type
diff_HC_MDD <- posterior_samples_latency[, "sub_type-HC"] - posterior_samples_latency[, "sub_type-MDD"]
# Calculate the Bayes factor for the hypothesis
bf_hypothesis <- mean(diff_HC_MDD > 0) / mean(diff_HC_MDD < 0)
# Print the hypothesis test results
print(bf_hypothesis) # 0.000, which suppport the hypothesis that HC in sub_type is different from MDD in sub_type
## [1] 2.891051
# Perform Bayesian ANOVA for emo_dwell
dwell_emo_data <- data2_emo_dwell_subed
# Ensure subnum is a factor
dwell_emo_data$subnum <- as.factor(dwell_emo_data$subnum)
dwell_emo_data$sub_type <- as.factor(dwell_emo_data$sub_type)
dwell_emo_data$expression <- as.factor(dwell_emo_data$expression)
# Perform Bayesian ANOVA
bf_anova_dwell_emo <- anovaBF(dwell_emo ~ sub_type * expression + subnum, data = dwell_emo_data, whichRandom = "subnum")
## Warning: data coerced from tibble to data frame
# Print the results
print(bf_anova_dwell_emo)
## Bayes factor analysis
## --------------
## [1] sub_type + subnum : 0.2978441 ±2.96%
## [2] expression + subnum : 0.1119341 ±0.57%
## [3] sub_type + expression + subnum : 0.03321896 ±1.59%
## [4] sub_type + expression + sub_type:expression + subnum : 0.0001727863 ±2.69%
##
## Against denominator:
## dwell_emo ~ subnum
## ---
## Bayes factor type: BFlinearModel, JZS
# Extract posterior samples from the full model (index 4)
# The full model includes sub_type, expression, their interaction, and subnum.
posterior_samples_dwell_emo <- posterior(bf_anova_dwell_emo, iterations = 1000, index = 4)
# summary(posterior_samples_dwell_emo)
# Test specific hypothesis: HC in sub_type is different from MDD in sub_type
# Calculate the difference between HC and MDD in sub_type
diff_HC_MDD <- posterior_samples_dwell_emo[, "sub_type-HC"] - posterior_samples_dwell_emo[, "sub_type-MDD"]
# Calculate the Bayes factor for the hypothesis
bf_hypothesis <- mean(diff_HC_MDD > 0) / mean(diff_HC_MDD < 0)
# Print the hypothesis test results
print(bf_hypothesis) # 0.000, which suppport the hypothesis that HC in sub_type is different from MDD in sub_type
## [1] 0.1520737
# Perform Bayesian ANOVA for neu_dwell
dwell_neu_data <- data2_neu_dwell_subed
# Ensure subnum is a factor
dwell_neu_data$subnum <- as.factor(dwell_neu_data$subnum)
dwell_neu_data$sub_type <- as.factor(dwell_neu_data$sub_type)
dwell_neu_data$expression <- as.factor(dwell_neu_data$expression)
# Perform Bayesian ANOVA
bf_anova_dwell_neu <- anovaBF(dwell_neu ~ sub_type * expression + subnum, data = dwell_neu_data, whichRandom = "subnum")
## Warning: data coerced from tibble to data frame
# Print the results
print(bf_anova_dwell_neu)
## Bayes factor analysis
## --------------
## [1] sub_type + subnum : 0.3632505 ±1.42%
## [2] expression + subnum : 0.06137201 ±1.54%
## [3] sub_type + expression + subnum : 0.02241055 ±1.56%
## [4] sub_type + expression + sub_type:expression + subnum : 0.0003384799 ±7.89%
##
## Against denominator:
## dwell_neu ~ subnum
## ---
## Bayes factor type: BFlinearModel, JZS
# Extract posterior samples from the full model (index 4)
# The full model includes sub_type, expression, their interaction, and subnum.
posterior_samples_dwell_neu <- posterior(bf_anova_dwell_neu, iterations = 1000, index = 4)
# summary(posterior_samples_dwell_neu)
# Test specific hypothesis: HC in sub_type is different from MDD in sub_type
# Calculate the difference between HC and MDD in sub_type
diff_HC_MDD <- posterior_samples_dwell_neu[, "sub_type-HC"] - posterior_samples_dwell_neu[, "sub_type-MDD"]
# Calculate the Bayes factor for the hypothesis
bf_hypothesis <- mean(diff_HC_MDD > 0) / mean(diff_HC_MDD < 0)
# Print the hypothesis test results
print(bf_hypothesis) # 0.1, which suppport the hypothesis that HC in sub_type is not different from MDD in sub_type
## [1] 0.08932462