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
sumrepdat3E_sub <- summarySE(data3, measurevar = "latency",
groupvars=c('subnum',"sub_type", "trial_type","expression",'AOI_emo')) %>%
mutate(IQR_L=outlier.IQR(latency_mean)[1], IQR_H=outlier.IQR(latency_mean)[2]) %>%
mutate(valid1=if_else(latency_mean>=IQR_L & latency_mean<=IQR_H, 1, 0)) %>%
filter(valid1==1) # need to be tested
## Warning in qt(conf.interval/2 + 0.5, datac$N - 1): NaNs produced
# Perform ANOVA for first fixation latency
anova_fix_latency <- aov_ez(
data = sumrepdat3E_sub,
id = 'subnum',
dv = 'latency_mean',
between = 'sub_type',
within = 'expression',
observed = 'sub_type'
)
## Converting to factor: sub_type
## Warning: More than one observation per design cell, aggregating data using `fun_aggregate = mean`.
## To turn off this warning, pass `fun_aggregate = mean` explicitly.
## 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 886.625 1038.0000 NaN
## # 64 MD039_ MDD NaN 925.1667 836.75
## 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) 174436262 1 5745442 122 3704.0189 < 2e-16 ***
## sub_type 256021 1 5745442 122 5.4364 0.02136 *
## expression 19191 2 1582109 244 1.4799 0.22970
## sub_type:expression 2598 2 1582109 244 0.2003 0.81859
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## expression 0.96997 0.15812
## sub_type:expression 0.96997 0.15812
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## expression 0.97085 0.2301
## sub_type:expression 0.97085 0.8122
##
## HF eps Pr(>F[HF])
## expression 0.9863076 0.2298797
## sub_type:expression 0.9863076 0.8156391
# Perform post hoc test on 'sub_type' variable
emmeans_latency_result <- emmeans(anova_fix_latency, pairwise ~ sub_type, adjust = "bonferroni")
# Print the post hoc test results
print(emmeans_latency_result)
## $emmeans
## sub_type emmean SE df lower.CL upper.CL
## HC 669 17.5 122 634 704
## MDD 722 14.7 122 693 751
##
## Results are averaged over the levels of: expression
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## HC - MDD -53.3 22.9 122 -2.332 0.0214
##
## Results are averaged over the levels of: expression
# 引导条件下首滞留
data2 <- my_data %>%
filter(trial_type != 'centered', continue_gaze==1) %>%
mutate(FIX_latency=endT-Tstart) %>%
group_by(subnum, Ntrial, sub_type, trial_type, expression) %>%
summarise(latency=max(FIX_latency)) %>%
ungroup() %>%
filter(latency>=120) %>% # need papers to support
mutate(latency=if_else(latency>=2500, 2500, latency))
## `summarise()` has grouped output by 'subnum', 'Ntrial', 'sub_type',
## 'trial_type'. You can override using the `.groups` argument.
sumrepdat2E_sub <- summarySE(data2, measurevar = "latency",
groupvars=c('subnum',"sub_type", "trial_type", "expression")) %>%
mutate(IQR_L=outlier.IQR(latency_mean)[1], IQR_H=outlier.IQR(latency_mean)[2]) %>%
mutate(valid1=if_else(latency_mean>=IQR_L & latency_mean<=IQR_H, 1, 0)) %>%
filter(valid1==1) # need to be tested
# Perform ANOVA for first fixation latency
anova_dwell <- aov_ez(
data = sumrepdat2E_sub,
id = 'subnum',
dv = 'latency_mean',
between = 'sub_type',
within = 'expression',
observed = 'sub_type'
)
## Converting to factor: sub_type
## Warning: More than one observation per design cell, aggregating data using `fun_aggregate = mean`.
## To turn off this warning, pass `fun_aggregate = mean` explicitly.
## Warning: Missing values for 3 ID(s), which were removed before analysis:
## HC081_, HC107_, MD070_
## Below the first few rows (in wide format) of the removed cases with missing data.
## subnum sub_type F H S
## # 29 HC081_ HC 891.8571 NaN 739.5938
## # 42 HC107_ HC NaN NaN 1046.0625
## # 86 MD070_ MDD NaN 949.0625 867.0000
## Contrasts set to contr.sum for the following variables: sub_type
# Summarize the ANOVA results
summary(anova_dwell)
##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 83422827 1 9311717 118 1057.1512 < 2.2e-16 ***
## sub_type 152687 1 9311717 118 1.9349 0.1668
## expression 110830 2 1352069 236 9.6725 9.174e-05 ***
## sub_type:expression 622 2 1352069 236 0.0543 0.9471
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## expression 0.89703 0.0017344
## sub_type:expression 0.89703 0.0017344
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## expression 0.90664 0.0001674 ***
## sub_type:expression 0.90664 0.9339428
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## expression 0.9199872 0.0001535728
## sub_type:expression 0.9199872 0.9360288921
# Perform post hoc test on 'expression' variable
emmeans_result <- emmeans(anova_dwell, pairwise ~ expression, adjust = "bonferroni")
# Print the post hoc test results
print(emmeans_result)
## $emmeans
## expression emmean SE df lower.CL upper.CL
## F 510 16.8 118 477 544
## H 492 16.3 118 460 524
## S 467 15.1 118 437 497
##
## 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 18.4 11.42 118 1.613 0.3284
## F - S 43.5 8.97 118 4.855 <.0001
## H - S 25.1 9.25 118 2.717 0.0228
##
## Results are averaged over the levels of: sub_type
## P value adjustment: bonferroni method for 3 tests