Loaded Joint Simon Effect

Effects of working memory load on the joint Simon effect

set.seed(12345) # for reproducibility

if (!require("pacman", quietly = TRUE)) install.packages("pacman")
pacman::p_load(tidyverse, psych, papaja, knitr, cowplot) # data manipulation, visualization
pacman::p_load(afex, permuco, Superpower) # stats
pacman::p_load_gh("thomasp85/patchwork", "RLesur/klippy", "mitchelloharawild/icons")

options(dplyr.summarise.inform=FALSE)
options(knitr.kable.NA = '')
set_sum_contrasts() # see Singmann & Kellen (2020)
klippy::klippy()

nsims <- 1e4

1 Task Comparison

아래 모든 과제에서 참가자는 화살표의 방향을 무시한 채로 색깔에 대해서만 반응한다. 가령, (A) 사이먼 과제에서 참가자는 화살표가 초록색이면 왼쪽 버튼을, 빨간색이면 오른쪽 버튼을 누른다. (B) Go/no-go 과제의 참가자는 화살표가 빨간색일 때에만 반응하고 초록색일 때는 반응하지 않아야 한다. (C) 결합 사이먼 과제에서 한 참가자는 초록색 화살표에 대해 go/no-go 과제를 하고, 다른 참가자는 빨간색 화살표에 대해 go/no-go 과제를 한다.

ggdraw() + draw_image("figs/Fig1.jpg")

---

2 Pilot Experiment

Sebanz 등(2003)이 보고한 효과를 재현하였다. 참가자는 Solo 구획에서는 과제를 혼자 수행하였고(Fig 1B), Joint 구획에서는 다른 참가자(실험자)와 함께 과제를 수행하였다(Fig 1C).

2.1 Data

전체 시행 중 참가자가 반응해야 하는 시행을 추출하여 분석하였다.

E1 = read.csv('data/mergeJSEv1.csv', header = TRUE)
E1$sid <- factor(E1$sid)
E1$performer = factor(E1$performer, levels=c(1,2), labels=c("Solo","Joint"))
E1$compatibility = factor(E1$compatibility, levels=c(1,2), 
                          labels=c("Compatible","Incompatible"))
E1$rt <- E1$rt*1000 
str(E1)
## 'data.frame':    4800 obs. of  9 variables:
##  $ sid          : Factor w/ 10 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ performer    : Factor w/ 2 levels "Solo","Joint": 1 1 1 1 1 1 1 1 1 1 ...
##  $ compatibility: Factor w/ 2 levels "Compatible","Incompatible": 2 1 2 2 1 1 2 1 2 1 ...
##  $ turn         : int  1 0 0 0 1 0 1 0 1 1 ...
##  $ dir          : int  1 1 2 2 2 1 1 1 1 2 ...
##  $ col          : int  1 2 2 2 1 2 1 2 1 1 ...
##  $ resp         : int  0 0 0 0 1 0 1 0 1 1 ...
##  $ corr         : int  0 1 1 1 1 1 1 1 1 1 ...
##  $ rt           : num  7000 7000 7000 7000 491 ...
headTail(E1)
##       sid performer compatibility turn dir col resp corr     rt
## 1       1      Solo  Incompatible    1   1   1    0    0   7000
## 2       1      Solo    Compatible    0   1   2    0    1   7000
## 3       1      Solo  Incompatible    0   2   2    0    1   7000
## 4       1      Solo  Incompatible    0   2   2    0    1   7000
## ...  <NA>      <NA>          <NA>  ... ... ...  ...  ...    ...
## 4797   10     Joint    Compatible    0   1   2    0    1   7000
## 4798   10     Joint    Compatible    1   2   1    1    1 447.23
## 4799   10     Joint  Incompatible    1   1   1    1    1 369.61
## 4800   10     Joint    Compatible    1   2   1    1    1 483.01

table(E1$sid)
## 
##   1   2   3   4   5   6   7   8   9  10 
## 480 480 480 480 480 480 480 480 480 480
table(E1$corr, E1$sid)
##    
##       1   2   3   4   5   6   7   8   9  10
##   0   1   0   2   6   0   1   0   1   2   1
##   1 479 480 478 474 480 479 480 479 478 479
table(E1$turn, E1$sid) # 0-confederate turn, 1-subject turn
##    
##       1   2   3   4   5   6   7   8   9  10
##   0 240 240 240 240 240 240 240 240 240 240
##   1 240 240 240 240 240 240 240 240 240 240

S1 <- E1 %>% 
  filter(turn==1) %>% 
  select(sid, performer, compatibility, corr, rt)
  
colSums(table(S1$sid, S1$corr))
##    0    1 
##    9 2391
headTail(S1)
##       sid performer compatibility corr     rt
## 1       1      Solo  Incompatible    0   7000
## 2       1      Solo    Compatible    1 491.41
## 3       1      Solo  Incompatible    1 362.35
## 4       1      Solo  Incompatible    1 398.84
## ...  <NA>      <NA>          <NA>  ...    ...
## 2397   10     Joint  Incompatible    1 347.43
## 2398   10     Joint    Compatible    1 447.23
## 2399   10     Joint  Incompatible    1 369.61
## 2400   10     Joint    Compatible    1 483.01
str(S1)
## 'data.frame':    2400 obs. of  5 variables:
##  $ sid          : Factor w/ 10 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ performer    : Factor w/ 2 levels "Solo","Joint": 1 1 1 1 1 1 1 1 1 1 ...
##  $ compatibility: Factor w/ 2 levels "Compatible","Incompatible": 2 1 2 2 1 2 2 1 1 1 ...
##  $ corr         : int  0 1 1 1 1 1 1 1 1 1 ...
##  $ rt           : num  7000 491 362 399 451 ...

# 개별참가자의 표적 색깔
S2 <- E1 %>% 
  filter(turn==1) %>% 
  select(sid, performer, compatibility, col, resp, corr, rt)
table(S2$col, S2$sid) # 1,2,5,6,9,10 - red, 3,4,7,8 - green
##    
##       1   2   3   4   5   6   7   8   9  10
##   1 240 240   0   0 240 240   0   0 240 240
##   2   0   0 240 240   0   0 240 240   0   0

# check errors, anticipatory rts, etc.
( nCorr <- S1 %>% filter(corr==1) %>% nrow() ) # correct RT
## [1] 2391
( nIncorr <- S1 %>% filter(corr!=1) %>% nrow() ) # incorrect RT
## [1] 9
nIncorr*100/nrow(S1)
## [1] 0.375

( numAntici <- S1 %>% filter(corr==1 & rt < 200) %>% nrow() ) # anticipatory correct RT
## [1] 0

# correct RTs
cS1 <- S1 %>% filter(corr == 1) %>% 
  group_by(sid, performer, compatibility) %>% 
  nest() %>% 
  mutate(lbound = map(data, ~mean(.$rt)-3*sd(.$rt)),
         ubound = map(data, ~mean(.$rt)+3*sd(.$rt))) %>%
  unnest(c(lbound, ubound)) %>%
  unnest(data) %>%
  mutate(Outlier = (rt < lbound)|(rt > ubound)) %>%
  filter(Outlier == FALSE) %>%
  ungroup() %>%
  select(sid, performer, compatibility, rt)

# percetage trimmed
(nCorr-nrow(cS1))*100/nrow(S1)
## [1] 1.5

오류율은 0.375%에 불과하므로 조건별 비교를 생략한다.

2.2 Reaction Times

반응시간 중 200ms보다 짧은 기대반응(anticipatory responses)은 한 건도 없었다. 각 참가자의 조건별 반응시간들 중 평균으로부터 3SD를 벗어난 1.5%의 반응시간들을 분석에서 제외하였다.

2.2.1 Descriptive Stats

cS1slong <- cS1 %>% group_by(sid, performer, compatibility) %>% 
  summarise(rt = mean(rt)) %>% 
  ungroup()

# 4 cells
cS1slong %>% group_by(performer, compatibility) %>% 
  summarise(M = mean(rt), SD = sd(rt)) %>% 
  ungroup() %>% 
  kable(digits = 4, caption = "Descriptive statistics: Performer x Compatibility")

Descriptive statistics: Performer x Compatibility

performer	compatibility	M	SD
Solo	Compatible	359.6976	31.8009
Solo	Incompatible	360.7449	31.0670
Joint	Compatible	340.6182	34.6561
Joint	Incompatible	349.6920	35.0724

# 2 cells by performer
cS1 %>% group_by(sid, performer) %>% 
  summarise(rt = mean(rt)) %>% 
  ungroup() %>% 
  group_by(performer) %>% 
  summarise(M = mean(rt), SD = sd(rt)) %>% 
  ungroup() %>% 
  kable(digits = 4, caption = "Descriptive statistics: Performer")

Descriptive statistics: Performer

performer	M	SD
Solo	360.2414	31.2522
Joint	345.1581	34.6292

# 2 cells by compatibility
cS1 %>% group_by(sid, compatibility) %>% 
  summarise(rt = mean(rt)) %>% 
  ungroup() %>% 
  group_by(compatibility) %>% 
  summarise(M = mean(rt), SD = sd(rt)) %>% 
  ungroup() %>% 
  kable(digits = 4, caption = "Descriptive statistics: Compatibility")

Descriptive statistics: Compatibility

compatibility	M	SD
Compatible	350.2170	27.0015
Incompatible	355.2173	27.6542

2.2.2 Plots

# 2x2 factorial plot
cS1wsci <- cS1slong %>% 
  wsci(id = "sid", 
       factors = c("performer", "compatibility"), 
       dv = "rt") %>% 
  rename("wsci" = "rt")

cS1mean <- cS1slong %>% group_by(performer, compatibility) %>%
  summarise(rt = mean(rt)) %>% 
  ungroup()

cS1g <- merge(cS1mean, cS1wsci, by = c("performer", "compatibility"), all = TRUE)

cS1swide <- cS1slong %>%  # wide format, needed for geom_segments.
  pivot_wider(id_cols='sid', names_from=c('performer', 'compatibility'), values_from='rt' )

G1 <- ggplot() + 
  geom_bar(data=cS1g, aes(x=performer, y=rt, fill=compatibility),
           stat="identity", width=0.7, color="black", position=position_dodge(.8)) + 
  geom_linerange(data=cS1g, aes(x=performer, ymin=rt-wsci, ymax=rt+wsci, group=compatibility),
                 size=0.8, position=position_dodge(0.8)) +
  scale_fill_manual(values=c('gray100','gray30'), # 회색조 
                    labels=c("Compatible", "Incompatible")) +
  # scale_fill_manual(values=c('#0073C2FF','#EFC000FF'), # 파란-노랑
  #                   labels=c("Compatible", "Incompatible")) +
  geom_point(data=cS1slong, aes(x=performer, y=rt, group=compatibility),
             position=position_dodge(0.6), color="gray80", size=1.8) +
  geom_segment(data=cS1swide, aes(x=1-.15, y=Solo_Compatible, 
                                  xend=1+.15, yend=Solo_Incompatible), color="gray80") +
  geom_segment(data=cS1swide, aes(x=2-.15, y=Joint_Compatible, 
                                  xend=2+.15, yend=Joint_Incompatible), color="gray80") +
  labs(x = "Performer", y = "Reaction Time (ms)") +
  coord_cartesian(ylim = c(250, 450), clip = "on") +
  theme_bw(base_size = 18) +
  theme(legend.position="top",
        legend.spacing.x = unit(0.5, 'lines'),
        legend.title = element_blank(),
        legend.background = element_blank(),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank())


# compatibility effect
cS1CEswide <- cS1swide %>% 
  mutate(Solo = Solo_Incompatible - Solo_Compatible,
         Joint = Joint_Incompatible - Joint_Compatible) %>% 
  select(sid, Solo, Joint)

cS1CEslong <- cS1CEswide %>% 
  pivot_longer(cols = c('Solo', 'Joint'), 
               names_to = 'performer',
               values_to = 'ce')
cS1CEslong$performer <- factor(cS1CEslong$performer, levels=c("Solo", "Joint"))

cS1CEmean <- cS1CEslong %>% group_by(performer) %>%
  summarise(ce = mean(ce)) %>% 
  ungroup()

cS1CEwsci <- cS1CEslong %>% # wsci - Morey
  wsci(id = "sid", 
       factors = c("performer"), 
       dv = "ce") %>% 
  rename("wsci" = "ce")

cS1CEg <- merge(cS1CEmean, cS1CEwsci, by = c("performer"), all = TRUE)

G2 <- ggplot(cS1CEslong, aes(x=performer, y=ce)) +
  geom_hline(yintercept = 0) +
  geom_violin(width = 0.5, trim=TRUE) + 
  geom_point(color="gray80", size=1.8) +
  geom_segment(data=cS1CEswide, color="gray80", 
               aes(x=1, y=Solo, xend=2, yend=Joint)) +
  geom_pointrange(cS1CEg, inherit.aes=FALSE,
                  mapping=aes(x = performer, y=ce, 
                              ymin = ce - wsci, ymax = ce + wsci), 
                  colour="darkred", size = 1) +
  labs(x = "Performer", y = "Congruency Effect (ms)") +
  coord_cartesian(ylim = c(-10, 30), clip = "on") +
  scale_y_continuous(breaks=c(-10, 0, 10, 20, 30)) +
  theme_bw(base_size = 18) +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank())

# Multipanel plot
G1 + G2 + plot_layout(nrow = 1, widths = c(2, 1.3))

2.2.3 Normality Test

반응시간 평균의 분포가 정규분포를 이루는지 검증하였다.

# normality test
N1 <- ggpubr::ggdensity(cS1slong$rt,
                        main = "Density plot",
                        xlab = "Response Times (msec)") + 
    theme(axis.text.y = element_blank(),
                axis.ticks.y = element_blank())
N2 <- ggpubr::ggqqplot(cS1slong$rt, 
                       main = "Q-Q plot")
N1 + N2 + plot_layout(nrow = 1, widths = c(1, 1))

shapiro.test(cS1slong$rt) # 만족

    Shapiro-Wilk normality test

data:  cS1slong$rt
W = 0.96988, p-value = 0.3567

정규분포 가정에 위배되지 않으므로 ANOVA를 실시할 수 있다.

2.2.4 ANOVA

cS1.aov <- cS1slong %>% aov_ez(id = "sid", dv = "rt", within = c("performer", "compatibility"))
cS1.aov %>% anova(es = "pes") %>% 
  kable(digits = 4, caption = "2x2 ANOVA table")

2x2 ANOVA table

	num Df	den Df	MSE	F	pes	Pr(>F)
performer	1	9	1402.4343	1.6185	0.1524	0.2352
compatibility	1	9	33.6637	7.6074	0.4581	0.0222
performer:compatibility	1	9	24.0489	6.6972	0.4266	0.0293

# afex_plot(cS1.aov, x = "performer", trace = "compatibility") 

cS1slong %>% filter(performer == 'Solo') %>% 
  aov_ez(id = 'sid', dv = 'rt', within = 'compatibility') %>% 
  anova(es = "pes") %>% 
  kable(digits = 4, caption = "Congruency effect in the Solo condition")

Congruency effect in the Solo condition

	num Df	den Df	MSE	F	pes	Pr(>F)
compatibility	1	9	23.4708	0.2337	0.0253	0.6404

cS1slong %>% filter(performer == 'Joint') %>% 
  aov_ez(id = 'sid', dv = 'rt', within = 'compatibility') %>% 
  anova(es = "pes") %>% 
  kable(digits = 4, caption = "Congruency effect in the Joint condition")

Congruency effect in the Joint condition

	num Df	den Df	MSE	F	pes	Pr(>F)
compatibility	1	9	34.2418	12.0224	0.5719	0.0071

2.2.5 Permutation Test

본실험의 자료는 정규성 가정에 위배된 경우가 있었다. 따라서 치환검정을 모든 분석에 적용하기로 한다.

# two-way permutation anova
cS1.2wayperm <- aovperm(rt ~ performer * compatibility + Error(sid/(performer*compatibility)),
                        data = cS1slong, np = nsims)
summary(cS1.2wayperm) %>%
    kable(digits = 4, caption = "Nonparametric ANOVA")

Nonparametric ANOVA

	SSn	dfn	SSd	dfd	MSEn	MSEd	F	parametric P(>F)	resampled P(>F)
performer	2269.8938	1	12621.9084	9	2269.8938	1402.4343	1.6185	0.2352	0.2411
compatibility	256.0934	1	302.9733	9	256.0934	33.6637	7.6074	0.0222	0.0241
performer:compatibility	161.0594	1	216.4402	9	161.0594	24.0489	6.6972	0.0293	0.0313

plot(cS1.2wayperm, effect = c("performer", 
                                                            "compatibility", 
                                                            "performer:compatibility"))

compatibility 주효과와 performer:compatibility 상호작용이 유의미하였다.
조건별 일치효과가 0보다 큰지 확인하였다.

cS1Solo.1wayperm <- aovperm(rt ~ compatibility + Error(sid/(compatibility)),
                        data = filter(cS1slong, performer=='Solo'), np = nsims)
summary(cS1Solo.1wayperm) %>%
    kable(digits = 4, caption = "Nonparametric simple test")

Nonparametric simple test

	SSn	dfn	SSd	dfd	MSEn	MSEd	F	parametric P(>F)	resampled P(>F)
compatibility	5.4847	1	211.2373	9	5.4847	23.4708	0.2337	0.6404	0.6386

cS1Joint.1wayperm <- aovperm(rt ~ compatibility + Error(sid/(compatibility)),
                            data = filter(cS1slong, performer=='Joint'), np = nsims)
summary(cS1Joint.1wayperm) %>%
    kable(digits = 4, caption = "Nonparametric simple test")

Nonparametric simple test

	SSn	dfn	SSd	dfd	MSEn	MSEd	F	parametric P(>F)	resampled P(>F)
compatibility	411.6681	1	308.1762	9	411.6681	34.2418	12.0224	0.0071	0.0078

사후검증을 위해 FDR 절차로 p-value를 보정하였다.

tmp <- p.adjust(c(cS1Solo.1wayperm$table$`resampled P(>F)`, 
                  cS1Joint.1wayperm$table$`resampled P(>F)`), "fdr")
data.frame(Solo = tmp[1], Joint = tmp[2]) %>% 
  kable(align = 'c', digits = 4, caption = "FDR adjusted p-values")

FDR adjusted p-values

Solo	Joint
0.6386	0.0156

일치 효과는 Solo 조건에서 유의미하지 않고 Joint 조건에서만 유의미하였다.

2.3 Power Test

예비실험의 Joint 조건의 결과를 활용하여 참가자 수를 추정하였다. High load 집단의 dual 과제에서만 SE가 관찰되지 않는 이원상호작용을 기대하였다.

# 조건별 평균
diffM <- cS1 %>% 
  group_by(sid, performer, compatibility) %>% 
  summarise(M = mean(rt)) %>% 
  ungroup() %>% 
  group_by(performer, compatibility) %>% 
  summarise(RT = mean(M)) %>% 
  ungroup() %>% 
    filter(performer == "Joint") %>% 
    pivot_wider(names_from = compatibility, values_from = RT) %>% 
    mutate(RTdiff = Incompatible - Compatible)
as.numeric(diffM$RTdiff)
[1] 9.073787

# 평균 표준편차
pooledSD <- cS1 %>% 
  group_by(sid, performer, compatibility) %>% 
  summarise(M = mean(rt)) %>% 
  ungroup() %>% 
  group_by(performer, compatibility) %>% 
  summarise(SD = sd(M)) %>% 
  ungroup() %>% 
    filter(performer == "Joint") %>%
    summarise(pooledSD = mean(SD))
as.numeric(pooledSD)
[1] 34.86426

# 집단별 참가자내 변인의 상관
wcor <- cS1 %>% 
  group_by(sid, performer, compatibility) %>% 
  summarise(M = mean(rt)) %>% 
  ungroup() %>%
    filter(performer == "Joint") %>%
  pivot_wider(names_from = compatibility, values_from = M) %>% 
    summarise(wcor = cor(Compatible, Incompatible))
as.numeric(wcor)
[1] 0.9718997

dx0 <- ANOVA_design(
    design = "2b*2w", 
    n = 10,
    mu = c(9,9,9,0),
    sd = 35,
    r = 0.97,
    labelnames = c("GROUP", "Lowload", "Highload",
                                 "TASK", "Single", "Dual"),
    plot = FALSE
)

plot_power(dx0, 
           min_n = 10, max_n = 50, 
           desired_power = 95, 
                     plot = TRUE)

Achieved Power and Sample Size for ANOVA-level effects
    variable                     label  n achieved_power desired_power
1      GROUP Desired Power Not Reached 50           9.83            95
2       TASK    Desired Power Achieved 25          95.32            95
3 GROUP:TASK    Desired Power Achieved 25          95.32            95

이원 상호작용 관찰을 위해 각 집단에 최소 25명이 필요하다. 본실험 참가자 수는 8의 배수가 되어야 하므로 Low load 집단과 High load 집단에 각 32명을 모집하기로 결정하였다.

---

3 Main Experiment

작업기억 부하 효과를 관찰하였다. 한 세션에 두 명의 참가자가 함께 과제를 수행하였다(Fig 1C).

3.1 Working Memory Task

3.1.1 Descriptive Stats

조건별 작업기억 과제의 정확률은 다음과 같다.

WW <- read.csv('data/mergeJSEv5_WMtask.csv', header = TRUE)
headTail(WW)
##     sid group trial change resp corr   rt
## 1   101   Low     1      1    1    1  0.9
## 2   101   Low     2      1    1    1 0.62
## 3   101   Low     3      1    1    1 0.69
## 4   101   Low     4      2    2    1 1.01
## ... ...  <NA>   ...    ...  ...  ...  ...
## 765 232  High     9      2    2    1 0.71
## 766 232  High    10      1    1    1 0.76
## 767 232  High    11      1    1    1 0.68
## 768 232  High    12      1    1    1 0.42

table(WW$corr, WW$sid)
##    
##     101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119
##   0   1   0   0   0   0   0   0   1   0   4   1   0   3   1   0   0   1   0   0
##   1  11  12  12  12  12  12  12  11  12   8  11  12   9  11  12  12  11  12  12
##    
##     120 121 122 123 124 125 126 127 128 129 130 131 132 201 202 203 204 205 206
##   0   0   0   1   1   0   0   0   1   1   3   0   1   3   1   1   6   1   1   1
##   1  12  12  11  11  12  12  12  11  11   9  12  11   9  11  11   6  11  11  11
##    
##     207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225
##   0   3   0   0   1   0   0   0   0   1   1   0   1   1   4   1   1   3   0   0
##   1   9  12  12  11  12  12  12  12  11  11  12  11  11   8  11  11   9  12  12
##    
##     226 227 228 229 230 231 232
##   0   3   3   3   1   3   0   1
##   1   9   9   9  11   9  12  11
WW <- WW %>% mutate(group = factor(group, levels = c('Low', 'High')))
WW$rt <- WW$rt*1000

# 참가자 수
WW %>% group_by(group) %>% 
  summarise(count = n_distinct(sid)) %>% 
  ungroup()
## # A tibble: 2 × 2
##   group count
##   <fct> <int>
## 1 Low      32
## 2 High     32

# WM RT
rWW <- WW %>% filter(corr == 1) %>% 
  group_by(sid, group) %>% 
  nest() %>% 
  mutate(lbound = map(data, ~mean(.$rt)-3*sd(.$rt)), # lower/upper bound 계산
         ubound = map(data, ~mean(.$rt)+3*sd(.$rt))) %>% 
  unnest(c(lbound, ubound)) %>% 
  unnest(data) %>% 
  mutate(outlier = (rt < lbound)|(rt > ubound)) %>% 
  filter(outlier == FALSE) %>% 
  ungroup() %>% 
  select(sid, group, rt)

(nrow(WW)-nrow(rWW))*100/nrow(WW)
## [1] 8.463542

aWWsum <- WW %>% group_by(group, sid) %>%
  summarise(MN = mean(corr)*100) %>% 
  ungroup() %>% 
  group_by(group) %>%
  summarise(N = n(),
                        Mean = mean(MN), 
            SD = sd(MN)) %>% 
  ungroup()

aWWsum %>% 
  kable(digits = 4, caption = "Descriptive statistics: WM accuracy")

Descriptive statistics: WM accuracy

group	N	Mean	SD
Low	32	94.0104	9.0385
High	32	89.0625	12.0553

평균으로부터 3SD 떨어진 반응시간 8.4635%를 제거한 후, 집단별 작업기억 과제의 반응시간은 다음과 같다.

rWWsum <- rWW %>% group_by(group, sid) %>% 
    summarise(MN = mean(rt)) %>% 
    ungroup() %>% 
    group_by(group) %>% 
    summarise(N = n(),
                        Mean = mean(MN), 
                        SD = sd(MN)) %>% 
    ungroup() 

rWWsum %>% 
    kable(digits = 4, caption = "Descriptive statistics: WM RT")

Descriptive statistics: WM RT

group	N	Mean	SD
Low	32	956.7012	153.3393
High	32	1134.0039	207.9728

3.1.2 Plots

aWWslong <- WW %>% group_by(sid, group) %>% 
  summarise(acc = mean(corr)*100) %>% 
  ungroup()

WG1 <- ggplot(aWWslong, aes(x=group, y=acc)) +
  geom_violin(width = 0.5, trim = TRUE) +
  geom_dotplot(binaxis='y', stackdir='center', dotsize=.7, 
               color='gray80', fill='gray80', binwidth=1.5) +
  stat_summary(fun.data = mean_cl_normal, geom="pointrange", color="darkred", size=1) +
  labs(x="Group", y="WM Accuracy") +
  scale_x_discrete(labels=c("Low load", "High load")) +
  coord_cartesian(ylim=c(50, 100)) +
  theme_bw(base_size = 14) +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        plot.title = element_text(hjust = .5))

# WM RT
rWWslong <- rWW %>% group_by(sid, group) %>% 
  summarise(rt = mean(rt)) %>% 
  ungroup()

WG2 <- ggplot(rWWslong, aes(x=group, y=rt)) +
  geom_violin(width = 0.5, trim = TRUE) +
  geom_dotplot(binaxis='y', stackdir='center', dotsize=.7, 
               color='gray80', fill='gray80', binwidth=30) +
  stat_summary(fun.data = mean_cl_normal, geom="pointrange", color="darkred", size=1) +
  labs(x="Group", y="Reaction Time (ms)") +
  scale_x_discrete(labels=c("Low load", "High load")) +
  coord_cartesian(ylim=c(500, 1500)) +
  theme_bw(base_size = 14) +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        plot.title = element_text(hjust = .5))

# Multipanel plot
WG1 + WG2

3.1.3 Accuracy Results

3.1.3.1 Normality Test

N3 <- ggpubr::ggdensity(aWWslong$acc, 
                        main = "Density plot",
                        xlab = "WM Accuracy") + 
    theme(axis.text.y = element_blank(),
                axis.ticks.y = element_blank())
N4 <- ggpubr::ggqqplot(aWWslong$acc, 
                       main = "Q-Q plot")
N3 + N4 + plot_layout(nrow = 1, widths = c(1, 1))

shapiro.test(aWWslong$acc)  # 가정 위배 

    Shapiro-Wilk normality test

data:  aWWslong$acc
W = 0.7388, p-value = 2.426e-09

3.1.3.2 Permutation Test

aWW.perm <- aovperm(acc ~ group, data = aWWslong, np = nsims)

Warning in check_distribution(distribution = distribution, digits = 10, : the
distribution of group may be discrete.

summary(aWW.perm) %>%
    kable(digits = 4, caption = "Nonparametric ANOVA")

Nonparametric ANOVA

	SS	df	F	parametric P(>F)	resampled P(>F)
group	391.7101	1	3.4508	0.068	0.0855
Residuals	7037.7604	62

치환검정에서는 두 집단의 정확도 차이가 유위미하지 않았다.

3.1.3.3 Power Test

daW <- ANOVA_design(
    design = "2b", 
    n = aWWsum$N,
    mu = aWWsum$Mean, 
    sd = aWWsum$SD,
    labelnames = c("GROUP", "Low", "High"),
    plot = FALSE
)

aWpwr <- ANOVA_power(daW, nsims = nsims)

Power and Effect sizes for ANOVA tests
            power effect_size
anova_GROUP 44.43     0.06532

Power and Effect sizes for pairwise comparisons (t-tests)
                       power effect_size
p_GROUP_Low_GROUP_High 44.43     -0.4693

시뮬레이션 결과, 정확도 차이의 검증력은 44.4%였다. 효과 크기는 평균 \(\eta^2_p\) = 0.065 또는 Cohen’s \(d\) = -0.469.

3.1.4 RT Results

3.1.4.1 Normality Test

N5 <- ggpubr::ggdensity(rWWslong$rt, 
                  main = "Density plot",
                  xlab = "WM Response Times (msec)")
N6 <- ggpubr::ggqqplot(rWWslong$rt, 
                       main = "Q-Q plot")
N5 + N6 + plot_layout(nrow = 1, widths = c(1, 1))

shapiro.test(rWWslong$rt)  #  정규분포다. 

## 
##  Shapiro-Wilk normality test
## 
## data:  rWWslong$rt
## W = 0.95972, p-value = 0.03531

반응시간 분포는 정규성 가정에 위배되지 않았다.

3.1.4.2 ANOVA

rWW.aov <- rWWslong %>% aov_ez(id = "sid", dv = "rt", between = "group")
rWW.aov %>% anova(es = "pes") %>% 
  kable(digits = 4, caption = "One-way ANOVA table")

One-way ANOVA table

	num Df	den Df	MSE	F	pes	Pr(>F)
group	1	62	33382.8	15.067	0.1955	3e-04

변량분석 결과, Low load 집단의 반응이 High load 집단보다 유의미하게 빨랐다.

3.1.4.3 Permutation Test

rWW.perm <- aovperm(rt ~ group, data = rWWslong, np = nsims)
summary(rWW.perm) %>%
    kable(digits = 4, caption = "Nonparametric ANOVA")

Nonparametric ANOVA

	SS	df	F	parametric P(>F)	resampled P(>F)
group	502980.1	1	15.067	3e-04	5e-04
Residuals	2069733.8	62

치환검정은 변량분석과 같은 결과를 보였다.

3.1.4.4 Power Test

drW <- ANOVA_design(
    design = "2b", 
    n = rWWsum$N,
    mu = rWWsum$Mean, 
    sd = rWWsum$SD,
    labelnames = c("GROUP", "Low", "High"),
    plot = FALSE
)

rWpwr <- ANOVA_power(drW, nsims = nsims)

Power and Effect sizes for ANOVA tests
            power effect_size
anova_GROUP 96.98      0.2024

Power and Effect sizes for pairwise comparisons (t-tests)
                       power effect_size
p_GROUP_Low_GROUP_High 96.98      0.9839

3.2 Color Judgment Task

본 연구에서 가장 중요한 분석이다.

참가자간 요인 Group (Low load vs. High load), 참가자내 요인 Task (Single vs. Dual), 참가자내 요인 Compatibility (Compatible vs. Incompatible)의 2x2x2 설계이다.

3.2.1 Descriptive Stats

3.2.1.1 Accuracy

TT <- read.csv('data/mergeJSEv5_GNGtask.csv', header = TRUE)
headTail(TT)
##       sid group initblk epoch trial task compatibility targ col dir resp corr
## 1     101   Low       1     1     1    1             1    1   1   1    1    1
## 2     101   Low       1     1     2    1             0    0   2   2    0    1
## 3     101   Low       1     1     3    1             0    0   2   1    0    1
## 4     101   Low       1     1     4    1             0    0   2   1    0    1
## ...   ...  <NA>     ...   ...   ...  ...           ...  ... ... ...  ...  ...
## 30717 232  High       1     4   477    2             0    0   1   2    0    1
## 30718 232  High       1     4   478    2             0    0   1   1    0    1
## 30719 232  High       1     4   479    2             1    1   2   2    1    1
## 30720 232  High       1     4   480    2             0    0   1   1    0    1
##         rt
## 1     0.99
## 2        0
## 3        0
## 4        0
## ...    ...
## 30717    0
## 30718    0
## 30719  0.3
## 30720    0

unique(TT$group)
## [1] "Low"  "High"
unique(TT$epoch)
## [1] 1 2 3 4
unique(TT$task)
## [1] 1 2
unique(TT$compatibility) # 1 0 2
## [1] 1 0 2
unique(TT$targ)
## [1] 1 0
unique(TT$corr)
## [1] 1 0

table(TT$targ, TT$sid)
##    
##     101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119
##   0 240 240 240 240 240 240 240 240 240 240 240 240 240 240 240 240 240 240 240
##   1 240 240 240 240 240 240 240 240 240 240 240 240 240 240 240 240 240 240 240
##    
##     120 121 122 123 124 125 126 127 128 129 130 131 132 201 202 203 204 205 206
##   0 240 240 240 240 240 240 240 240 240 240 240 240 240 240 240 240 240 240 240
##   1 240 240 240 240 240 240 240 240 240 240 240 240 240 240 240 240 240 240 240
##    
##     207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225
##   0 240 240 240 240 240 240 240 240 240 240 240 240 240 240 240 240 240 240 240
##   1 240 240 240 240 240 240 240 240 240 240 240 240 240 240 240 240 240 240 240
##    
##     226 227 228 229 230 231 232
##   0 240 240 240 240 240 240 240
##   1 240 240 240 240 240 240 240
table(TT$compatibility, TT$sid)
##    
##     101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119
##   0 240 240 240 240 240 240 240 240 240 240 240 240 240 240 240 240 240 240 240
##   1 120 120 120 120 120 120 120 120 120 120 120 120 120 120 120 120 120 120 120
##   2 120 120 120 120 120 120 120 120 120 120 120 120 120 120 120 120 120 120 120
##    
##     120 121 122 123 124 125 126 127 128 129 130 131 132 201 202 203 204 205 206
##   0 240 240 240 240 240 240 240 240 240 240 240 240 240 240 240 240 240 240 240
##   1 120 120 120 120 120 120 120 120 120 120 120 120 120 120 120 120 120 120 120
##   2 120 120 120 120 120 120 120 120 120 120 120 120 120 120 120 120 120 120 120
##    
##     207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225
##   0 240 240 240 240 240 240 240 240 240 240 240 240 240 240 240 240 240 240 240
##   1 120 120 120 120 120 120 120 120 120 120 120 120 120 120 120 120 120 120 120
##   2 120 120 120 120 120 120 120 120 120 120 120 120 120 120 120 120 120 120 120
##    
##     226 227 228 229 230 231 232
##   0 240 240 240 240 240 240 240
##   1 120 120 120 120 120 120 120
##   2 120 120 120 120 120 120 120


TT <- TT %>% mutate(group = factor(group, levels = c('Low', 'High')),
                    task = factor(task, levels=1:2, labels=c("Single","Dual")),
                    compatibility = factor(compatibility, levels=0:2, 
                                           labels=c("Nontarget","Compatible","Incompatible")),
                    targ = factor(targ, levels=0:1, labels=c("Nontarget","Target")),
                    rt = rt*1000) %>% 
  filter(targ == "Target") %>% 
  droplevels() %>% 
  select(sid, group, epoch, task, compatibility, corr, rt)
  
TT %>% sapply(levels)
## $sid
## NULL
## 
## $group
## [1] "Low"  "High"
## 
## $epoch
## NULL
## 
## $task
## [1] "Single" "Dual"  
## 
## $compatibility
## [1] "Compatible"   "Incompatible"
## 
## $corr
## NULL
## 
## $rt
## NULL

headTail(TT)
##       sid group epoch   task compatibility corr     rt
## 1     101   Low     1 Single    Compatible    1 988.24
## 2     101   Low     1 Single  Incompatible    1 326.42
## 3     101   Low     1 Single    Compatible    1 311.76
## 4     101   Low     1 Single  Incompatible    1 268.35
## ...   ...  <NA>   ...   <NA>          <NA>  ...    ...
## 15357 232  High     4   Dual    Compatible    1 269.62
## 15358 232  High     4   Dual  Incompatible    1 289.87
## 15359 232  High     4   Dual    Compatible    1  274.2
## 15360 232  High     4   Dual    Compatible    1 299.84

unique(TT$compatibility)
## [1] Compatible   Incompatible
## Levels: Compatible Incompatible

# averaged accuracy
aTTslong <- TT %>% group_by(group, sid, task, compatibility) %>% 
  summarise(Accuracy = mean(corr)*100) %>% 
  ungroup() 

# summary
aTTsum <- aTTslong %>% group_by(group, task, compatibility) %>% 
  summarise(N = n(),
                    MN = mean(Accuracy), 
            SD = sd(Accuracy)) %>% 
  ungroup() 
aTTsum %>% 
  kable(digits = 4, caption = "Descriptive statistics: Group x Load x Congruency")

Descriptive statistics: Group x Load x Congruency

group	task	compatibility	N	MN	SD
Low	Single	Compatible	32	97.1354	4.9047
Low	Single	Incompatible	32	96.7708	4.8810
Low	Dual	Compatible	32	97.6042	3.8085
Low	Dual	Incompatible	32	96.7188	4.7798
High	Single	Compatible	32	98.2292	3.6890
High	Single	Incompatible	32	97.0833	3.9937
High	Dual	Compatible	32	97.1354	5.9904
High	Dual	Incompatible	32	97.4479	5.2017

aTTslong %>% group_by(group) %>% 
  summarise(MN = mean(Accuracy), 
            SD = sd(Accuracy)) %>% 
  ungroup() %>% 
  kable(digits = 4, caption = "Descriptive statistics: Group")

Descriptive statistics: Group

group	MN	SD
Low	97.0573	4.5750
High	97.4740	4.7739

aTTslong %>% group_by(task) %>% 
  summarise(MN = mean(Accuracy), 
            SD = sd(Accuracy)) %>% 
  ungroup() %>% 
  kable(digits = 4, caption = "Descriptive statistics: Load")

Descriptive statistics: Load

task	MN	SD
Single	97.3047	4.3828
Dual	97.2266	4.9595

aTTslong %>% group_by(compatibility) %>% 
  summarise(MN = mean(Accuracy), 
            SD = sd(Accuracy)) %>% 
  ungroup() %>% 
  kable(digits = 4, caption = "Descriptive statistics: Congruency")

Descriptive statistics: Congruency

compatibility	MN	SD
Compatible	97.5260	4.6580
Incompatible	97.0052	4.6878

3.2.1.2 RT

( num_anticip <- nrow(TT %>% filter(corr==1 & rt < 150)) ) # anticipatory response 제거
## [1] 4

cTT <- TT %>% filter(corr==1 & rt >= 150) # 정반응만 선별

num_anticip*100/nrow(cTT) # 0.02678093%
## [1] 0.02678093

headTail(cTT)
##       sid group epoch   task compatibility corr     rt
## 1     101   Low     1 Single    Compatible    1 988.24
## 2     101   Low     1 Single  Incompatible    1 326.42
## 3     101   Low     1 Single    Compatible    1 311.76
## 4     101   Low     1 Single  Incompatible    1 268.35
## ...   ...  <NA>   ...   <NA>          <NA>  ...    ...
## 14933 232  High     4   Dual    Compatible    1 269.62
## 14934 232  High     4   Dual  Incompatible    1 289.87
## 14935 232  High     4   Dual    Compatible    1  274.2
## 14936 232  High     4   Dual    Compatible    1 299.84

# trimmed
tTT <- cTT %>% group_by(group, sid, task, compatibility) %>% 
  nest() %>% 
  mutate(lbound = map(data, ~mean(.$rt)-3*sd(.$rt)), # lower/upper bound 계산
         ubound = map(data, ~mean(.$rt)+3*sd(.$rt))) %>% 
  unnest(c(lbound, ubound)) %>% 
  unnest(data) %>% 
  mutate(outlier = (rt < lbound)|(rt > ubound)) %>% 
  filter(outlier == FALSE) %>% 
  select(group, sid, epoch, task, compatibility, rt)

(nrow(cTT)-nrow(tTT))*100/nrow(cTT) # 1.131494% 제거되었다.
## [1] 1.131494

tTTslong <- tTT %>% group_by(group, sid, task, compatibility) %>% 
  summarise(RT = mean(rt)) %>% 
  ungroup()

# summary
tTTsum <- tTTslong %>% group_by(group, task, compatibility) %>% 
  summarise(MN = mean(RT), 
            SD = sd(RT)) %>% 
  ungroup() 
tTTsum %>% 
  kable(digits = 4, caption = "Descriptive statistics: Group x Task x Compatibility")

Descriptive statistics: Group x Task x Compatibility

group	task	compatibility	MN	SD
Low	Single	Compatible	355.8836	26.2698
Low	Single	Incompatible	362.0058	29.0895
Low	Dual	Compatible	361.8038	26.5217
Low	Dual	Incompatible	366.3566	27.9480
High	Single	Compatible	354.2354	33.5891
High	Single	Incompatible	354.9457	34.3319
High	Dual	Compatible	359.3876	38.6664
High	Dual	Incompatible	360.8769	40.1932

tTTslong %>% group_by(group) %>% 
  summarise(MN = mean(RT), 
            SD = sd(RT)) %>% 
  ungroup() %>% 
  kable(digits = 4, caption = "Descriptive statistics: Group")

Descriptive statistics: Group

group	MN	SD
Low	361.5125	27.4105
High	357.3614	36.4755

tTTslong %>% group_by(task) %>% 
  summarise(MN = mean(RT), 
            SD = sd(RT)) %>% 
  ungroup() %>% 
  kable(digits = 4, caption = "Descriptive statistics: Task")

Descriptive statistics: Task

task	MN	SD
Single	356.7677	30.7842
Dual	362.1062	33.5924

tTTslong %>% group_by(compatibility) %>% 
  summarise(MN = mean(RT), 
            SD = sd(RT)) %>% 
  ungroup() %>% 
  kable(digits = 4, caption = "Descriptive statistics: Compatibility")

Descriptive statistics: Compatibility

compatibility	MN	SD
Compatible	357.8276	31.4529
Incompatible	361.0463	33.1055

반응시간이 150ms보다 빠른 4시행( 0.0268%)을 분석에서 제외한 후, 각 참가자의 조건별 반응시간들 중 평균으로부터 3SD를 벗어난 1.1315%의 반응시간들을 분석에서 추가로 제외하였다.

3.2.2 Plot

3.2.2.1 Accuracy

aTTswide <- aTTslong %>% pivot_wider(id_cols=c('group', 'sid'), 
                                     names_from=c('task', 'compatibility'), 
                                     values_from='Accuracy' )
aTTswideL <- aTTswide %>% filter(group == 'Low')
aTTswideH <- aTTswide %>% filter(group == 'High')

aTTmean <- aTTslong %>% 
  group_by(group, task, compatibility) %>% 
  summarise(Accuracy = mean(Accuracy)) %>% 
  ungroup() 

tmp1 <- aTTslong %>% filter(group == "Low") %>% 
  wsci(id = "sid",
       factor = c("task", "compatibility"),
       dv = "Accuracy") %>% 
  mutate(group = "Low") %>% 
  select(group, task, compatibility, Accuracy) %>% 
  rename("wsci" = "Accuracy")

tmp2 <- aTTslong %>% filter(group == "High") %>% 
  wsci(id = "sid",
       factor = c("task", "compatibility"),
       dv = "Accuracy") %>% 
  mutate(group = "High") %>% 
  select(group, task, compatibility, Accuracy) %>% 
  rename("wsci" = "Accuracy")

aTTwsci <- merge(tmp1, tmp2, all = TRUE)

aTTg <- merge(aTTmean, aTTwsci, by = c("group", "task", "compatibility"), all = TRUE)

group.labs <- c("Low load group", "High load group")
names(group.labs) <- c("Low", "High")

ggplot() + 
  geom_bar(data=aTTg, aes(x=task, y=Accuracy, fill=compatibility),
           stat="identity", width=0.7, color="black", position=position_dodge(.8)) + 
  facet_wrap(~group, labeller = labeller(group = group.labs)) +
  geom_linerange(data=aTTg, aes(x=task, ymin=Accuracy-wsci, ymax=Accuracy+wsci, 
                                group=compatibility),
                 size=1, position=position_dodge(0.8)) +
  scale_fill_manual(values=c('gray100','gray30'),
                    labels=c("Compatible", "Incompatible")) +
  # scale_fill_manual(values=c('#0073C2FF','#EFC000FF'),
  #                   labels=c("Compatible", "Incompatible")) +
  geom_point(data=aTTslong, aes(x=task, y=Accuracy, group=compatibility),
             position=position_dodge(0.6), color="gray80", size=1.8) +
  geom_segment(data=aTTswideL, aes(x=1-.15, y=Single_Compatible,
                                   xend=1+.15, yend=Single_Incompatible), color="gray80") +
  geom_segment(data=aTTswideL, aes(x=2-.15, y=Dual_Compatible,
                                   xend=2+.15, yend=Dual_Incompatible), color="gray80") +
  geom_segment(data=aTTswideH, aes(x=1-.15, y=Single_Compatible,
                                   xend=1+.15, yend=Single_Incompatible), color="gray80") +
  geom_segment(data=aTTswideH, aes(x=2-.15, y=Dual_Compatible,
                                   xend=2+.15, yend=Dual_Incompatible), color="gray80") +
  labs(x = "Task", y = "Accuracy") +
  coord_cartesian(ylim = c(50, 100), clip = "on") +
  theme_bw(base_size = 18) +
  theme(legend.position="top",
        legend.spacing.x = unit(0.5, 'lines'),
        strip.text.x = element_text(size = 18),
        legend.title = element_blank(),
        legend.background = element_blank(),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank())

3.2.2.2 RT

tTTswide <- tTTslong %>% pivot_wider(id_cols=c('group', 'sid'), 
                                     names_from=c('task', 'compatibility'), 
                                     values_from='RT' )
tTTswideL <- tTTswide %>% filter(group == 'Low')
tTTswideH <- tTTswide %>% filter(group == 'High')

tTTmean <- tTTslong %>% 
  group_by(group, task, compatibility) %>% 
  summarise(RT = mean(RT)) %>% 
  ungroup() 

tmp1 <- tTTslong %>% filter(group == "Low") %>% 
  wsci(id = "sid",
       factor = c("task", "compatibility"),
       dv = "RT") %>% 
  mutate(group = "Low") %>% 
  select(group, task, compatibility, RT) %>% 
  rename("wsci" = "RT")

tmp2 <- tTTslong %>% filter(group == "High") %>% 
  wsci(id = "sid",
       factor = c("task", "compatibility"),
       dv = "RT") %>% 
  mutate(group = "High") %>% 
  select(group, task, compatibility, RT) %>% 
  rename("wsci" = "RT")

tTTwsci <- merge(tmp1, tmp2, all = TRUE)

tTTg <- merge(tTTmean, tTTwsci, by = c("group", "task", "compatibility"), all = TRUE)

group.labs <- c("Low load group", "High load group")
names(group.labs) <- c("Low", "High")

# range(tTTslong$RT)

ggplot() + 
  geom_bar(data=tTTg, aes(x=task, y=RT, fill=compatibility),
           stat="identity", width=0.7, color="black", position=position_dodge(.8)) + 
  facet_wrap(~group, labeller = labeller(group = group.labs)) +
  geom_linerange(data=tTTg, aes(x=task, ymin=RT-wsci, ymax=RT+wsci, group=compatibility),
                 size=1, position=position_dodge(0.8)) +
  scale_fill_manual(values=c('gray100','gray30'),
                    labels=c("Compatible", "Incompatible")) +
  # scale_fill_manual(values=c('#0073C2FF','#EFC000FF'),
  #                   labels=c("Compatible", "Incompatible")) +
  geom_point(data=tTTslong, aes(x=task, y=RT, group=compatibility),
             position=position_dodge(0.6), color="gray80", size=1.8) +
  geom_segment(data=tTTswideL, aes(x=1-.15, y=Single_Compatible,
                                   xend=1+.15, yend=Single_Incompatible), color="gray80") +
  geom_segment(data=tTTswideL, aes(x=2-.15, y=Dual_Compatible,
                                   xend=2+.15, yend=Dual_Incompatible), color="gray80") +
  geom_segment(data=tTTswideH, aes(x=1-.15, y=Single_Compatible,
                                   xend=1+.15, yend=Single_Incompatible), color="gray80") +
  geom_segment(data=tTTswideH, aes(x=2-.15, y=Dual_Compatible,
                                   xend=2+.15, yend=Dual_Incompatible), color="gray80") +
  labs(x = "Task", y = "Reaction Times (msec)") +
  coord_cartesian(ylim = c(300, 550), clip = "on") +
  theme_bw(base_size = 18) +
  theme(legend.position="top",
        legend.spacing.x = unit(0.5, 'lines'),
        strip.text.x = element_text(size = 18),
        legend.title = element_blank(),
        legend.background = element_blank(),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank())

tTTCEswide <- tTTswide %>% mutate(absCE = Single_Incompatible - Single_Compatible,
                                  preCE = Dual_Incompatible - Dual_Compatible) %>% 
  select(group, sid, absCE, preCE)

tTTCEslong <- tTTCEswide %>% pivot_longer(cols = c('absCE', 'preCE'), names_to = "task")

tTTCEswideL <- tTTCEswide %>% filter(group == 'Low')
tTTCEswideH <- tTTCEswide %>% filter(group == 'High')

tTTCEmean <- tTTCEslong %>% 
  group_by(group, task) %>% 
  summarise(value = mean(value)) %>% 
  ungroup() 

tmp1 <- tTTCEslong %>% filter(group == "Low") %>% 
  wsci(id = "sid",
       factor = "task",
       dv = "value") %>% 
  mutate(group = "Low") %>% 
  select(group, task, value) %>% 
  rename("wsci" = "value")

tmp2 <- tTTCEslong %>% filter(group == "High") %>% 
  wsci(id = "sid",
       factor = "task",
       dv = "value") %>% 
  mutate(group = "High") %>% 
  select(group, task, value) %>% 
  rename("wsci" = "value")

tTTCEwsci <- merge(tmp1, tmp2, all = TRUE)

tTTCEg <- merge(tTTCEmean, tTTCEwsci, by = c("group", "task"), all = TRUE)

ggplot(data=tTTCEslong, aes(x=group, y=value, color=task)) +
  geom_hline(yintercept = 0) +
  geom_violin(width = 0.5, size=1, trim=TRUE) +
  geom_point(aes(x=group, y=value, group=task),
             position=position_dodge(0.5), color="gray80", size=1.8, show.legend = FALSE) +
  geom_segment(data=filter(tTTCEswide, group=="Low"), inherit.aes = FALSE,
               aes(x=1-.12, y=filter(tTTCEswide, group=="Low")$absCE,
                   xend=1+.12, yend=filter(tTTCEswide, group=="Low")$preCE),
               color="gray80") +
  geom_segment(data=filter(tTTCEswide, group=="High"), inherit.aes = FALSE,
               aes(x=2-.12, y=filter(tTTCEswide, group=="High")$absCE,
                   xend=2+.12, yend=filter(tTTCEswide, group=="High")$preCE),
               color="gray80") +
  geom_pointrange(data=tTTCEg, 
                  aes(x = group, ymin = value-wsci, ymax = value+wsci, group = task),
                  position = position_dodge(0.5), color = "darkred", size = 1, show.legend = FALSE) +
  scale_color_manual(values=c('#0073C2FF','#EFC000FF'),
                    labels=c("Single", "Dual")) +
  scale_x_discrete(labels=c("Low" = "Low load", "High" = "High load")) +
  labs(x = "Group", 
       y = "Compatibility Effect (ms)", 
       color='Task') +
  coord_cartesian(ylim = c(-20, 40), clip = "on") +
  theme_bw(base_size = 18) +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) 

3.2.3 Accuracy Results

3.2.3.1 Noramlity Test

N7 <- ggpubr::ggdensity(aTTslong$Accuracy, 
                        main = "Density plot",
                        xlab = "WM Accuracy")
N8 <- ggpubr::ggqqplot(aTTslong$Accuracy, 
                       main = "Q-Q plot")
N7 + N8 + plot_layout(nrow = 1, widths = c(1, 1))

shapiro.test(aTTslong$Accuracy)  # 가정 위배 

    Shapiro-Wilk normality test

data:  aTTslong$Accuracy
W = 0.64282, p-value < 2.2e-16

자료가 정규성 가정에 위배된다.

3.2.3.2 Permutation Test

aTT.perm <- aovperm(Accuracy ~ group * task * compatibility + Error(sid/(task*compatibility)),
                    data = aTTslong, np = nsims)
summary(aTT.perm) %>%
    kable(digits = 4, caption = "Nonparametric ANOVA")

Nonparametric ANOVA

	SSn	dfn	SSd	dfd	MSEn	MSEd	F	parametric P(>F)	resampled P(>F)
group	11.1111	1	4641.4931	62	11.1111	74.8628	0.1484	0.7014	0.7083
task	0.3906	1	391.5799	62	0.3906	6.3158	0.0618	0.8044	0.8041
group:task	5.2517	1	391.5799	62	5.2517	6.3158	0.8315	0.3654	0.3638
compatibility	17.3611	1	234.7222	62	17.3611	3.7858	4.5858	0.0362	0.0344
group:compatibility	0.6944	1	234.7222	62	0.6944	3.7858	0.1834	0.6699	0.6731
task:compatibility	3.5156	1	241.9271	62	3.5156	3.9020	0.9010	0.3462	0.3559
group:task:compatibility	15.6684	1	241.9271	62	15.6684	3.9020	4.0154	0.0495	0.0525

compatibility 주효과만 유의미하였다. group:load:cong 상호작용이 유의미한 경향성을 보였다.

aTTLS.perm <- aovperm(Accuracy ~ compatibility + Error(sid/(compatibility)),
                      data = filter(aTTslong, group=='Low', task=='Single'), np = nsims)
summary(aTTLS.perm) %>%
    kable(digits = 4, caption = "Low load & Single task: Simon effect")

Low load & Single task: Simon effect

	SSn	dfn	SSd	dfd	MSEn	MSEd	F	parametric P(>F)	resampled P(>F)
compatibility	2.1267	1	143.7066	31	2.1267	4.6357	0.4588	0.5032	0.5197

aTTLD.perm <- aovperm(Accuracy ~ compatibility + Error(sid/(compatibility)),
                      data = filter(aTTslong, group=='Low', task=='Dual'), np = nsims)
summary(aTTLD.perm) %>%
    kable(digits = 4, caption = "Low load & Dual task: Simon effect")

Low load & Dual task: Simon effect

	SSn	dfn	SSd	dfd	MSEn	MSEd	F	parametric P(>F)	resampled P(>F)
compatibility	12.5434	1	105.5122	31	12.5434	3.4036	3.6853	0.0641	0.0681

aTTHS.perm <- aovperm(Accuracy ~ compatibility + Error(sid/(compatibility)),
                      data = filter(aTTslong, group=='High', task=='Single'), np = nsims)
summary(aTTHS.perm) %>%
    kable(digits = 4, caption = "High load & Single task: Simon effect")

High load & Single task: Simon effect

	SSn	dfn	SSd	dfd	MSEn	MSEd	F	parametric P(>F)	resampled P(>F)
compatibility	21.0069	1	81.7708	31	21.0069	2.6378	7.9639	0.0083	0.0068

aTTHD.perm <- aovperm(Accuracy ~ compatibility + Error(sid/(compatibility)),
                      data = filter(aTTslong, group=='High', task=='Dual'), np = nsims)
summary(aTTHD.perm) %>%
    kable(digits = 4, caption = "High load & Dual task: Simon effect")

High load & Dual task: Simon effect

	SSn	dfn	SSd	dfd	MSEn	MSEd	F	parametric P(>F)	resampled P(>F)
compatibility	1.5625	1	145.6597	31	1.5625	4.6987	0.3325	0.5683	0.5707

사후검증에서 FDR 절차로 p-value를 보정하였다.

tmp <- p.adjust(c(aTTLS.perm$table$`resampled P(>F)`, 
                  aTTLD.perm$table$`resampled P(>F)`,
                  aTTHS.perm$table$`resampled P(>F)`,
                  aTTHD.perm$table$`resampled P(>F)`), "fdr")
data.frame(LowloadSolo = tmp[1], LowloadDual = tmp[2],
           HighloadSolo = tmp[3], HighloadDual = tmp[4]) %>% 
  kable(align = 'c', digits = 4, caption = "FDR adjusted p-values")

FDR adjusted p-values

LowloadSolo	LowloadDual	HighloadSolo	HighloadDual
0.5707	0.1362	0.0272	0.5707

3.2.3.3 Power Test

daT <- ANOVA_design(
    design = "2b*2w*2w", 
    n = 32,
    mu = aTTsum$MN, 
    sd = aTTsum$SD,
    labelnames = c("G", "Lo", "Hi",
                                 "T", "Sing", "Dual",
                                 "C", "Cmp", "Inc"),
    plot = FALSE
)

daTpwr <- ANOVA_power(daT, verbose = FALSE, nsims = nsims)

daTpwr$main_results %>%
    kable(digits = 4, caption = "Effect Sizes: ANOVA")

Effect Sizes: ANOVA

	power	effect_size
anova_G	10.22	0.0231
anova_T	5.11	0.0161
anova_G:T	7.85	0.0196
anova_C	14.48	0.0279
anova_G:C	5.08	0.0162
anova_T:C	7.16	0.0184
anova_G:T:C	13.56	0.0266

daTpwr$pc_results[c(1,14,23,28),] %>%
    kable(digits = 4, caption = "Effect Sizes: Post-hoc")

Effect Sizes: Post-hoc

	power	effect_size
p_G_Lo_T_Sing_C_Cmp_G_Lo_T_Sing_C_Inc	5.94	-0.0758
p_G_Lo_T_Dual_C_Cmp_G_Lo_T_Dual_C_Inc	12.38	-0.2083
p_G_Hi_T_Sing_C_Cmp_G_Hi_T_Sing_C_Inc	21.89	-0.3043
p_G_Hi_T_Dual_C_Cmp_G_Hi_T_Dual_C_Inc	5.80	0.0555

3.2.4 RT Results

3.2.4.1 Normality Test

N9 <- ggpubr::ggdensity(tTTslong$RT, 
                        main = "Density plot",
                        xlab = "Attention Task RT")
N10 <- ggpubr::ggqqplot(tTTslong$RT, 
                       main = "Q-Q plot")
N9 + N10 + plot_layout(nrow = 1, widths = c(1, 1))

shapiro.test(tTTslong$RT) # 가정 위배 

    Shapiro-Wilk normality test

data:  tTTslong$RT
W = 0.91493, p-value = 6.715e-11

반응시간 자료도 정규성 가정에 위배되었다.

3.2.4.2 Permutation Test

tTT.perm <- aovperm(RT ~ group * task * compatibility + Error(sid/(task*compatibility)),
                    data = tTTslong, np = nsims)
summary(tTT.perm) %>%
    kable(digits = 4, caption = "Nonparametric ANOVA")

Nonparametric ANOVA

	SSn	dfn	SSd	dfd	MSEn	MSEd	F	parametric P(>F)	resampled P(>F)
group	1102.8047	1	245617.772	62	1102.8047	3961.5770	0.2784	0.5997	0.6060
task	1824.0043	1	9822.264	62	1824.0043	158.4236	11.5135	0.0012	0.0010
group:task	2.6397	1	9822.264	62	2.6397	158.4236	0.0167	0.8977	0.9028
compatibility	663.0186	1	3790.207	62	663.0186	61.1324	10.8456	0.0016	0.0018
group:compatibility	287.3328	1	3790.207	62	287.3328	61.1324	4.7002	0.0340	0.0328
task:compatibility	2.4994	1	2356.338	62	2.4994	38.0055	0.0658	0.7985	0.7973
group:task:compatibility	22.0592	1	2356.338	62	22.0592	38.0055	0.5804	0.4490	0.4545

plot(tTT.perm, effect = c("group", "task", "compatibility"))

plot(tTT.perm, effect = c("group:task", "group:compatibility"))

plot(tTT.perm, effect = c("task:compatibility", "group:task:compatibility") )

조건별 일치효과가 0보다 큰지 확인하였다.

tTTLS.perm <- aovperm(RT ~ compatibility + Error(sid/(compatibility)),
                      data = filter(tTTslong, group=='Low', task=='Single'), np = nsims)
summary(tTTLS.perm) %>%
    kable(digits = 4, caption = "Low load & Single task: Simon effect")

Low load & Single task: Simon effect

	SSn	dfn	SSd	dfd	MSEn	MSEd	F	parametric P(>F)	resampled P(>F)
compatibility	599.7039	1	1903.501	31	599.7039	61.4033	9.7666	0.0038	0.0043

tTTLD.perm <- aovperm(RT ~ compatibility + Error(sid/(compatibility)),
                      data = filter(tTTslong, group=='Low', task=='Dual'), np = nsims)
summary(tTTLD.perm) %>%
    kable(digits = 4, caption = "Low load & Dual task: Simon effect")

Low load & Dual task: Simon effect

	SSn	dfn	SSd	dfd	MSEn	MSEd	F	parametric P(>F)	resampled P(>F)
compatibility	331.6474	1	1670.003	31	331.6474	53.8711	6.1563	0.0187	0.0209

tTTHS.perm <- aovperm(RT ~ compatibility + Error(sid/(compatibility)),
                      data = filter(tTTslong, group=='High', task=='Single'), np = nsims)
summary(tTTHS.perm) %>%
    kable(digits = 4, caption = "High load & Single task: Simon effect")

High load & Single task: Simon effect

	SSn	dfn	SSd	dfd	MSEn	MSEd	F	parametric P(>F)	resampled P(>F)
compatibility	8.0726	1	1191.88	31	8.0726	38.4478	0.21	0.65	0.6449

tTTHD.perm <- aovperm(RT ~ compatibility + Error(sid/(compatibility)),
                      data = filter(tTTslong, group=='High', task=='Dual'), np = nsims)
summary(tTTHD.perm) %>%
    kable(digits = 4, caption = "High load & Dual task: Simon effect")

High load & Dual task: Simon effect

	SSn	dfn	SSd	dfd	MSEn	MSEd	F	parametric P(>F)	resampled P(>F)
compatibility	35.486	1	1381.161	31	35.486	44.5536	0.7965	0.379	0.392

사후검증에서 FDR 절차로 p-value를 보정하였다.

tmp <- p.adjust(c(tTTLS.perm$table$`resampled P(>F)`, 
                  tTTLD.perm$table$`resampled P(>F)`,
                  tTTHS.perm$table$`resampled P(>F)`,
                  tTTHD.perm$table$`resampled P(>F)`), "fdr")
data.frame(LowloadSolo = tmp[1], LowloadDual = tmp[2],
           HighloadSolo = tmp[3], HighloadDual = tmp[4]) %>% 
  kable(align = 'c', digits = 4, caption = "FDR adjusted p-values")

FDR adjusted p-values

LowloadSolo	LowloadDual	HighloadSolo	HighloadDual
0.0172	0.0418	0.6449	0.5227

3.2.4.3 Power Test

drT <- ANOVA_design(
    design = "2b*2w*2w", 
    n = 32,
    mu = tTTsum$MN, 
    sd = tTTsum$SD,
    labelnames = c("G", "Lo", "Hi",
                                 "T", "Sing", "Dual",
                                 "C", "Cmp", "Inc"),
    plot = FALSE
)

drTpwr <- ANOVA_power(drT, verbose = FALSE, nsims = nsims)

drTpwr$main_results %>%
    kable(digits = 4, caption = "Effect Sizes: ANOVA")

Effect Sizes: ANOVA

	power	effect_size
anova_G	17.79	0.0319
anova_T	25.18	0.0412
anova_G:T	4.98	0.0160
anova_C	12.46	0.0255
anova_G:C	8.71	0.0204
anova_T:C	5.02	0.0159
anova_G:T:C	4.84	0.0160

drTpwr$pc_results[c(1,14,23,28),] %>%
    kable(digits = 4, caption = "Effect Sizes: Post-hoc")

Effect Sizes: Post-hoc

	power	effect_size
p_G_Lo_T_Sing_C_Cmp_G_Lo_T_Sing_C_Inc	14.06	0.2231
p_G_Lo_T_Dual_C_Cmp_G_Lo_T_Dual_C_Inc	10.61	0.1719
p_G_Hi_T_Sing_C_Cmp_G_Hi_T_Sing_C_Inc	5.17	0.0225
p_G_Hi_T_Dual_C_Cmp_G_Hi_T_Dual_C_Inc	5.41	0.0364

---

4 Plots

4.1 Experiment 1

G3 <- ggplot() + 
  geom_bar(data=cS1g, aes(x=performer, y=rt, fill=compatibility),
           stat="identity", width=0.7, color="black", position=position_dodge(.8)) + 
  geom_linerange(data=cS1g, aes(x=performer, ymin=rt-wsci, ymax=rt+wsci, group=compatibility),
                 size=0.8, position=position_dodge(0.8)) +
  scale_fill_manual(values=c('gray100','gray30'),
                    labels=c("Compatible", "Incompatible")) +
  labs(x = "Performer", y = "Reaction Time (ms)") +
  coord_cartesian(ylim = c(300, 400), clip = "on") +
  theme_bw(base_size = 16) +
  theme(legend.position="top",
        legend.spacing.x = unit(0.5, 'lines'),
        legend.title = element_blank(),
        legend.background = element_blank(),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank())

G4 <- ggplot(cS1CEslong, aes(x=performer, y=ce)) +
  geom_hline(yintercept = 0) +
  geom_violin(width = 0.5, trim=TRUE) + 
  ggbeeswarm::geom_quasirandom(color = "black", size = 3, alpha = 0.2, width = 0.2) +
  geom_pointrange(cS1CEg, inherit.aes=FALSE,
                  mapping=aes(x = performer, y=ce, 
                              ymin = ce - wsci, ymax = ce + wsci), 
                  colour="darkred", size = 1) +
  labs(x = "Performer", y = "Incompatible - Compatible") +
  coord_cartesian(ylim = c(-10, 30), clip = "on") +
  scale_y_continuous(breaks=c(-10, 0, 10, 20, 30)) +
  theme_bw(base_size = 16) +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank())

# Multipanel plot
G3 + G4 + plot_layout(nrow = 1, widths = c(2, 1.3))

4.2 Experiment 2

ggplot() + 
  geom_bar(data=tTTg, aes(x=task, y=RT, fill=compatibility),
           stat="identity", width=0.7, color="black", position=position_dodge(.8)) + 
  facet_wrap(~group, labeller = labeller(group = group.labs)) +
  geom_linerange(data=tTTg, aes(x=task, ymin=RT-wsci, ymax=RT+wsci, group=compatibility),
                 size=1, position=position_dodge(0.8)) +
  scale_fill_manual(values=c('gray100','gray30'),
                    labels=c("Compatible", "Incompatible")) +
  labs(x = "Task", y = "Reaction Times (msec)") +
  coord_cartesian(ylim = c(300, 400), clip = "on") +
  theme_bw(base_size = 16) +
  theme(legend.position="top",
        legend.spacing.x = unit(0.5, 'lines'),
        strip.text.x = element_text(size = 18),
        legend.title = element_blank(),
        legend.background = element_blank(),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank())

ggplot(data=tTTCEslong, aes(x=group, y=value, color=task)) +
  geom_hline(yintercept = 0) +
  geom_violin(width = 0.5, size=1, trim=TRUE) +
  ggbeeswarm::geom_quasirandom(groupOnX = TRUE, dodge.width=.5, 
                               size = 3, alpha = 0.2, width = 0.2, 
                               show.legend = FALSE) +
  geom_pointrange(data=tTTCEg, 
                  aes(x = group, ymin = value-wsci, ymax = value+wsci, group = task),
                  position = position_dodge(0.5), color = "darkred", size = 1, show.legend = FALSE) +
  scale_color_manual(values=c('#0073C2FF','#EFC000FF'),
                     labels=c("Single", "Dual")) +
  scale_x_discrete(labels=c("Low" = "Low load", "High" = "High load")) +
  labs(x = "Group", 
       y = "Incompatible - Compatible", 
       color='Task') +
  coord_cartesian(ylim = c(-20, 40), clip = "on") +
  theme_bw(base_size = 16) +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank())

---

5 Session Info

sessionInfo()
## R version 4.1.3 (2022-03-10)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur/Monterey 10.16
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] icons_0.2.0       klippy_0.0.0.9500 patchwork_1.1.1   Superpower_0.1.2 
##  [5] permuco_1.1.1     afex_1.1-0        lme4_1.1-29       Matrix_1.4-1     
##  [9] cowplot_1.1.1     knitr_1.38        papaja_0.1.0.9997 psych_2.2.3      
## [13] forcats_0.5.1     stringr_1.4.0     dplyr_1.0.8       purrr_0.3.4      
## [17] readr_2.1.2       tidyr_1.2.0       tibble_3.1.6      ggplot2_3.3.5    
## [21] tidyverse_1.3.1   pacman_0.5.1     
## 
## loaded via a namespace (and not attached):
##   [1] ggbeeswarm_0.6.0    TH.data_1.1-0       minqa_1.2.4        
##   [4] colorspace_2.0-3    ggsignif_0.6.3      ellipsis_0.3.2     
##   [7] htmlTable_2.4.0     estimability_1.3    base64enc_0.1-3    
##  [10] fs_1.5.2            rstudioapi_0.13     ggpubr_0.4.0       
##  [13] farver_2.1.0        fansi_1.0.3         mvtnorm_1.1-3      
##  [16] lubridate_1.8.0     xml2_1.3.3          codetools_0.2-18   
##  [19] splines_4.1.3       mnormt_2.0.2        Formula_1.2-4      
##  [22] jsonlite_1.8.0      nloptr_2.0.0        broom_0.8.0        
##  [25] cluster_2.1.3       dbplyr_2.1.1        png_0.1-7          
##  [28] compiler_4.1.3      httr_1.4.2          emmeans_1.7.3      
##  [31] backports_1.4.1     assertthat_0.2.1    fastmap_1.1.0      
##  [34] cli_3.2.0           htmltools_0.5.2     tools_4.1.3        
##  [37] lmerTest_3.1-3      coda_0.19-4         gtable_0.3.0       
##  [40] glue_1.6.2          reshape2_1.4.4      rappdirs_0.3.3     
##  [43] Rcpp_1.0.8.3        carData_3.0-5       cellranger_1.1.0   
##  [46] jquerylib_0.1.4     vctrs_0.4.1         nlme_3.1-157       
##  [49] xfun_0.30           rvest_1.0.2         lifecycle_1.0.1    
##  [52] rstatix_0.7.0       MASS_7.3-56         zoo_1.8-9          
##  [55] scales_1.2.0        hms_1.1.1           parallel_4.1.3     
##  [58] sandwich_3.0-1      RColorBrewer_1.1-3  yaml_2.3.5         
##  [61] gridExtra_2.3       sass_0.4.1          rpart_4.1.16       
##  [64] latticeExtra_0.6-29 stringi_1.7.6       highr_0.9          
##  [67] checkmate_2.0.0     permute_0.9-7       boot_1.3-28        
##  [70] rlang_1.0.2         pkgconfig_2.0.3     evaluate_0.15      
##  [73] lattice_0.20-45     labeling_0.4.2      htmlwidgets_1.5.4  
##  [76] tidyselect_1.1.2    plyr_1.8.7          magrittr_2.0.3     
##  [79] R6_2.5.1            magick_2.7.3        generics_0.1.2     
##  [82] Hmisc_4.6-0         multcomp_1.4-18     DBI_1.1.2          
##  [85] foreign_0.8-82      pillar_1.7.0        haven_2.4.3        
##  [88] withr_2.5.0         nnet_7.3-17         survival_3.3-1     
##  [91] abind_1.4-5         modelr_0.1.8        crayon_1.5.1       
##  [94] car_3.0-12          utf8_1.2.2          tmvnsim_1.0-2      
##  [97] tzdb_0.3.0          rmarkdown_2.13      jpeg_0.1-9         
## [100] grid_4.1.3          readxl_1.4.0        data.table_1.14.2  
## [103] reprex_2.0.1        digest_0.6.29       xtable_1.8-4       
## [106] numDeriv_2016.8-1.1 munsell_0.5.0       beeswarm_0.4.0     
## [109] vipor_0.4.5         bslib_0.3.1

---

View on Github