Define sample

Note: The analyses use the sample of subset subjects.

For all subjects, 8 pilot subjects were included. These were tested before formal preregistration.
For subset subjects, 8 pilot subjects were excluded.

tmpaEEG=read.csv(file.path(dir_res,"data_audio_signal_noise_N_33.tsv"), sep="\t",header=TRUE, dec=".")
#head(tmpaEEG)
# Analysis of visual P3
tmpvEEG=read.csv(file.path(dir_res,"data_visualP3_N_33.tsv"), sep="\t",header=TRUE, dec=".")
#head(tmpwmc)
tmpbeh=read.csv(file.path(dir_res,"data_logfiles.tsv"), sep="\t",header=TRUE, dec=".")
#head(tmpbeh)
tmpnasa=read.csv(file.path(dir_res,"data_NASA.tsv"), sep="\t",header=TRUE, dec=".")
#head(tmpnasa)
D=tmpbeh
D=merge(D, tmpaEEG, by = "fp")
D=merge(D, tmpnasa, by = "fp")
D=merge(D, tmpvEEG, by = "fp", all.x = TRUE)
rm(list=ls(pattern="tmp"))
if (any(as.numeric(table(D$fp))!=1)){
  stop("Subject numbers are not unique!")
  }
if (nowsample == "subset") {
  excludefps = 5:12
  indx = D$fp %in% excludefps
  D = D[!indx,]
  rm(indx)}
  #table(D$fp)

Quick analyses

Sample size: N = 25.

Sanity check for behavioral performance (across low and high load, and across frequencies):
Min number of hits (could be 0) : 48
Max number of hits (could be 96) : 96
Min number of false alarms (could be 0) : 0
Max number of false alarms (could be 302) : 68
Note that hits are correct responses between 200 and 1000 ms after target onset (defined in ASSR2_process_beh.R)
False alarms are responses at other times: 199*2 (number of nontarget trials times 2 blocks) - 96 (trial after a target) = 302.

D$HRlo20Hzcorr = (D$nHRlo20Hz + 0.5)/(96+1)
D$HRlo40Hzcorr = (D$nHRlo40Hz + 0.5)/(96+1)
D$HRlo80Hzcorr = (D$nHRlo80Hz + 0.5)/(96+1)
D$HRhi20Hzcorr = (D$nHRhi20Hz + 0.5)/(96+1)
D$HRhi40Hzcorr = (D$nHRhi40Hz + 0.5)/(96+1)
D$HRhi80Hzcorr = (D$nHRhi80Hz + 0.5)/(96+1)
# 48*2 = 96 targets
D$FAlo20Hzcorr = (D$nFAlo20Hz + 0.5)/(398+1-96)
D$FAlo40Hzcorr = (D$nFAlo40Hz + 0.5)/(398+1-96)
D$FAlo80Hzcorr = (D$nFAlo80Hz + 0.5)/(398+1-96)
D$FAhi20Hzcorr = (D$nFAhi20Hz + 0.5)/(398+1-96)
D$FAhi40Hzcorr = (D$nFAhi40Hz + 0.5)/(398+1-96)
D$FAhi80Hzcorr = (D$nFAhi80Hz + 0.5)/(398+1-96)
# 199*2 = 398 nontargets
# do not count the next trial after a target (because this is considered to be a valid response)
# Dprime
D$dprlo20Hz = qnorm(D$HRlo20Hzcorr) - qnorm(D$FAlo20Hzcorr)
D$dprlo40Hz = qnorm(D$HRlo40Hzcorr) - qnorm(D$FAlo40Hzcorr)
D$dprlo80Hz = qnorm(D$HRlo80Hzcorr) - qnorm(D$FAlo80Hzcorr)
D$dprhi20Hz = qnorm(D$HRhi20Hzcorr) - qnorm(D$FAhi20Hzcorr)
D$dprhi40Hz = qnorm(D$HRhi40Hzcorr) - qnorm(D$FAhi40Hzcorr)
D$dprhi80Hz = qnorm(D$HRhi80Hzcorr) - qnorm(D$FAhi80Hzcorr)
# do not compute criterion
#tmp$criterion = - (1/2) * (qnorm(tmp$aware_crit_corrected) + qnorm(tmp$aware_catch_corrected))
D$'dprlo-hi20Hz' = D$dprlo20Hz - D$dprhi20Hz
D$'dprlo-hi40Hz' = D$dprlo40Hz - D$dprhi40Hz
D$'dprlo-hi80Hz' = D$dprlo80Hz - D$dprhi80Hz
D$'meanRTmslo-hi20Hz' = D$meanRTmslo20Hz - D$meanRTmshi20Hz
D$'meanRTmslo-hi40Hz' = D$meanRTmslo40Hz - D$meanRTmshi40Hz
D$'meanRTmslo-hi80Hz' = D$meanRTmslo80Hz - D$meanRTmshi80Hz
# the numbers 1 to 2 refer to different blocks

# Amp S-N raw
D$AmpSmNno20Hz = rowMeans(cbind(D$AmpSmNno20_1,D$AmpSmNno20_2))
D$AmpSmNno40Hz = rowMeans(cbind(D$AmpSmNno40_1,D$AmpSmNno40_2))
D$AmpSmNno80Hz = rowMeans(cbind(D$AmpSmNno80_1,D$AmpSmNno80_2))
D$AmpSmNlo20Hz = rowMeans(cbind(D$AmpSmNlo20_1,D$AmpSmNlo20_2))
D$AmpSmNlo40Hz = rowMeans(cbind(D$AmpSmNlo40_1,D$AmpSmNlo40_2))
D$AmpSmNlo80Hz = rowMeans(cbind(D$AmpSmNlo80_1,D$AmpSmNlo80_2))
D$AmpSmNhi20Hz = rowMeans(cbind(D$AmpSmNhi20_1,D$AmpSmNhi20_2))
D$AmpSmNhi40Hz = rowMeans(cbind(D$AmpSmNhi40_1,D$AmpSmNhi40_2))
D$AmpSmNhi80Hz = rowMeans(cbind(D$AmpSmNhi80_1,D$AmpSmNhi80_2))

# Amp S raw
D$AmpSno20Hz = rowMeans(cbind(D$AmpSno20_1,D$AmpSno20_2))
D$AmpSno40Hz = rowMeans(cbind(D$AmpSno40_1,D$AmpSno40_2))
D$AmpSno80Hz = rowMeans(cbind(D$AmpSno80_1,D$AmpSno80_2))
D$AmpSlo20Hz = rowMeans(cbind(D$AmpSlo20_1,D$AmpSlo20_2))
D$AmpSlo40Hz = rowMeans(cbind(D$AmpSlo40_1,D$AmpSlo40_2))
D$AmpSlo80Hz = rowMeans(cbind(D$AmpSlo80_1,D$AmpSlo80_2))
D$AmpShi20Hz = rowMeans(cbind(D$AmpShi20_1,D$AmpShi20_2))
D$AmpShi40Hz = rowMeans(cbind(D$AmpShi40_1,D$AmpShi40_2))
D$AmpShi80Hz = rowMeans(cbind(D$AmpShi80_1,D$AmpShi80_2))

# Amp N raw
D$AmpNno20Hz = rowMeans(cbind(D$AmpNno20_1,D$AmpNno20_2))
D$AmpNno40Hz = rowMeans(cbind(D$AmpNno40_1,D$AmpNno40_2))
D$AmpNno80Hz = rowMeans(cbind(D$AmpNno80_1,D$AmpNno80_2))
D$AmpNlo20Hz = rowMeans(cbind(D$AmpNlo20_1,D$AmpNlo20_2))
D$AmpNlo40Hz = rowMeans(cbind(D$AmpNlo40_1,D$AmpNlo40_2))
D$AmpNlo80Hz = rowMeans(cbind(D$AmpNlo80_1,D$AmpNlo80_2))
D$AmpNhi20Hz = rowMeans(cbind(D$AmpNhi20_1,D$AmpNhi20_2))
D$AmpNhi40Hz = rowMeans(cbind(D$AmpNhi40_1,D$AmpNhi40_2))
D$AmpNhi80Hz = rowMeans(cbind(D$AmpNhi80_1,D$AmpNhi80_2))

# Amp SNR raw
D$AmpSNRno20Hz = rowMeans(cbind(D$AmpSNRno20_1,D$AmpSNRno20_2))
D$AmpSNRno40Hz = rowMeans(cbind(D$AmpSNRno40_1,D$AmpSNRno40_2))
D$AmpSNRno80Hz = rowMeans(cbind(D$AmpSNRno80_1,D$AmpSNRno80_2))
D$AmpSNRlo20Hz = rowMeans(cbind(D$AmpSNRlo20_1,D$AmpSNRlo20_2))
D$AmpSNRlo40Hz = rowMeans(cbind(D$AmpSNRlo40_1,D$AmpSNRlo40_2))
D$AmpSNRlo80Hz = rowMeans(cbind(D$AmpSNRlo80_1,D$AmpSNRlo80_2))
D$AmpSNRhi20Hz = rowMeans(cbind(D$AmpSNRhi20_1,D$AmpSNRhi20_2))
D$AmpSNRhi40Hz = rowMeans(cbind(D$AmpSNRhi40_1,D$AmpSNRhi40_2))
D$AmpSNRhi80Hz = rowMeans(cbind(D$AmpSNRhi80_1,D$AmpSNRhi80_2))
# the numbers 1 to 2 refer to different blocks

# ITC S-N raw
D$ItcSmNno20Hz = rowMeans(cbind(D$ItcSmNno20_1,D$ItcSmNno20_2))
D$ItcSmNno40Hz = rowMeans(cbind(D$ItcSmNno40_1,D$ItcSmNno40_2))
D$ItcSmNno80Hz = rowMeans(cbind(D$ItcSmNno80_1,D$ItcSmNno80_2))
D$ItcSmNlo20Hz = rowMeans(cbind(D$ItcSmNlo20_1,D$ItcSmNlo20_2))
D$ItcSmNlo40Hz = rowMeans(cbind(D$ItcSmNlo40_1,D$ItcSmNlo40_2))
D$ItcSmNlo80Hz = rowMeans(cbind(D$ItcSmNlo80_1,D$ItcSmNlo80_2))
D$ItcSmNhi20Hz = rowMeans(cbind(D$ItcSmNhi20_1,D$ItcSmNhi20_2))
D$ItcSmNhi40Hz = rowMeans(cbind(D$ItcSmNhi40_1,D$ItcSmNhi40_2))
D$ItcSmNhi80Hz = rowMeans(cbind(D$ItcSmNhi80_1,D$ItcSmNhi80_2))

# ITC S raw
D$ItcSno20Hz = rowMeans(cbind(D$ItcSno20_1,D$ItcSno20_2))
D$ItcSno40Hz = rowMeans(cbind(D$ItcSno40_1,D$ItcSno40_2))
D$ItcSno80Hz = rowMeans(cbind(D$ItcSno80_1,D$ItcSno80_2))
D$ItcSlo20Hz = rowMeans(cbind(D$ItcSlo20_1,D$ItcSlo20_2))
D$ItcSlo40Hz = rowMeans(cbind(D$ItcSlo40_1,D$ItcSlo40_2))
D$ItcSlo80Hz = rowMeans(cbind(D$ItcSlo80_1,D$ItcSlo80_2))
D$ItcShi20Hz = rowMeans(cbind(D$ItcShi20_1,D$ItcShi20_2))
D$ItcShi40Hz = rowMeans(cbind(D$ItcShi40_1,D$ItcShi40_2))
D$ItcShi80Hz = rowMeans(cbind(D$ItcShi80_1,D$ItcShi80_2))

# ITC N raw
D$ItcNno20Hz = rowMeans(cbind(D$ItcNno20_1,D$ItcNno20_2))
D$ItcNno40Hz = rowMeans(cbind(D$ItcNno40_1,D$ItcNno40_2))
D$ItcNno80Hz = rowMeans(cbind(D$ItcNno80_1,D$ItcNno80_2))
D$ItcNlo20Hz = rowMeans(cbind(D$ItcNlo20_1,D$ItcNlo20_2))
D$ItcNlo40Hz = rowMeans(cbind(D$ItcNlo40_1,D$ItcNlo40_2))
D$ItcNlo80Hz = rowMeans(cbind(D$ItcNlo80_1,D$ItcNlo80_2))
D$ItcNhi20Hz = rowMeans(cbind(D$ItcNhi20_1,D$ItcNhi20_2))
D$ItcNhi40Hz = rowMeans(cbind(D$ItcNhi40_1,D$ItcNhi40_2))
D$ItcNhi80Hz = rowMeans(cbind(D$ItcNhi80_1,D$ItcNhi80_2))

# ITC SNR raw
D$ItcSNRno20Hz = rowMeans(cbind(D$ItcSNRno20_1,D$ItcSNRno20_2))
D$ItcSNRno40Hz = rowMeans(cbind(D$ItcSNRno40_1,D$ItcSNRno40_2))
D$ItcSNRno80Hz = rowMeans(cbind(D$ItcSNRno80_1,D$ItcSNRno80_2))
D$ItcSNRlo20Hz = rowMeans(cbind(D$ItcSNRlo20_1,D$ItcSNRlo20_2))
D$ItcSNRlo40Hz = rowMeans(cbind(D$ItcSNRlo40_1,D$ItcSNRlo40_2))
D$ItcSNRlo80Hz = rowMeans(cbind(D$ItcSNRlo80_1,D$ItcSNRlo80_2))
D$ItcSNRhi20Hz = rowMeans(cbind(D$ItcSNRhi20_1,D$ItcSNRhi20_2))
D$ItcSNRhi40Hz = rowMeans(cbind(D$ItcSNRhi40_1,D$ItcSNRhi40_2))
D$ItcSNRhi80Hz = rowMeans(cbind(D$ItcSNRhi80_1,D$ItcSNRhi80_2))
# visual P3 is difference between targets and nontargets.
# the numbers 1 to 4 refer to different blocks

D$vP3lo20 = rowMeans(cbind(D$vP3_lo_tar20_1 - D$vP3_lo_non20_1, D$vP3_lo_tar20_2 - D$vP3_lo_non20_2))
D$vP3lo40 = rowMeans(cbind(D$vP3_lo_tar40_1 - D$vP3_lo_non40_1, D$vP3_lo_tar40_2 - D$vP3_lo_non40_2))
D$vP3lo80 = rowMeans(cbind(D$vP3_lo_tar80_1 - D$vP3_lo_non80_1, D$vP3_lo_tar80_2 - D$vP3_lo_non80_2))

D$vP3hi20 = rowMeans(cbind(D$vP3_hi_tar20_1 - D$vP3_hi_non20_1, D$vP3_hi_tar20_2 - D$vP3_hi_non20_2))
D$vP3hi40 = rowMeans(cbind(D$vP3_hi_tar40_1 - D$vP3_hi_non40_1, D$vP3_hi_tar40_2 - D$vP3_hi_non40_2))
D$vP3hi80 = rowMeans(cbind(D$vP3_hi_tar80_1 - D$vP3_hi_non80_1, D$vP3_hi_tar80_2 - D$vP3_hi_non80_2))

D$'vP3lo-hi20' = D$vP3lo20 - D$vP3hi20
D$'vP3lo-hi40' = D$vP3lo40 - D$vP3hi40
D$'vP3lo-hi80' = D$vP3lo80 - D$vP3hi80

Quick look at behavioral performance:
Simple hit rates (%) were as follows:
20-Hz AM:
Low load = 95.11
High load = 78.12
40-Hz AM:
Low load = 95.81
High load = 81.09
80-Hz AM:
Low load = 94.58
High load = 81.18

Manipulation checks of load

lo = low load
hi = high load
lo-hi = low minus high load

20, 40, 80 Hz = amplitude modulation frequency

LL and UL refer to the 95% confidence interval (from two-tailed t tests).

d prime and reaction time

dpr = d prime
meanRTms = RT to hits (ms)
Note: Because there was no task during no load, there were no targets and nontargets to examine d prime and RT.

table

Dtmp = subset(D, select=c('dprlo20Hz', 'dprhi20Hz', 'dprlo-hi20Hz',
                          'dprlo40Hz', 'dprhi40Hz', 'dprlo-hi40Hz',
                          'dprlo80Hz', 'dprhi80Hz', 'dprlo-hi80Hz',
                          'meanRTmslo20Hz', 'meanRTmshi20Hz', 'meanRTmslo-hi20Hz',
                          'meanRTmslo40Hz', 'meanRTmshi40Hz', 'meanRTmslo-hi40Hz',
                          'meanRTmslo80Hz', 'meanRTmshi80Hz', 'meanRTmslo-hi80Hz'))
RmCI = tableCIs(Dtmp)
RmCI %>%
  kable(digits = 3) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), 
                full_width = F, position = 'left')
Variable Mean LL UL N
dprlo20Hz 4.390 4.017 4.764 25
dprhi20Hz 2.557 2.253 2.860 25
dprlo-hi20Hz 1.834 1.484 2.183 25
dprlo40Hz 4.325 4.046 4.603 25
dprhi40Hz 2.629 2.344 2.914 25
dprlo-hi40Hz 1.696 1.445 1.947 25
dprlo80Hz 4.226 3.917 4.534 25
dprhi80Hz 2.558 2.296 2.820 25
dprlo-hi80Hz 1.668 1.383 1.953 25
meanRTmslo20Hz 372.097 349.208 394.985 25
meanRTmshi20Hz 486.641 462.119 511.162 25
meanRTmslo-hi20Hz -114.544 -137.326 -91.761 25
meanRTmslo40Hz 374.745 350.767 398.722 25
meanRTmshi40Hz 496.929 469.977 523.881 25
meanRTmslo-hi40Hz -122.184 -143.314 -101.055 25
meanRTmslo80Hz 370.060 352.275 387.844 25
meanRTmshi80Hz 494.460 470.406 518.513 25
meanRTmslo-hi80Hz -124.400 -144.095 -104.704 25
#(RmCI, file = paste0(file.path(dir_res,"table_descr_performance_",nowsample,"_subjects.tsv")), row.names = F, dec = ".", sep = "\t")

figure

Dtmp = subset(D, select=c('dprlo-hi20Hz', 'dprlo-hi40Hz', 'dprlo-hi80Hz'))
myylim = c(min(Dtmp), max(Dtmp)) 
RmCI = tableCIs(Dtmp)
#RmCI
Ddpr = RmCI[1:3,]
Dtmp = cbind(D$fp, Dtmp)
Dtmp = na.omit(Dtmp)
colnames(Dtmp) = c('fp','lo-hi20Hz','lo-hi40Hz','lo-hi80Hz')
bp1 = plotoverlay(dm = Dtmp, lbl = "d´", dvy = "Sensitivity (d´)", myylim)

Dtmp = subset(D, select=c('meanRTmslo-hi20Hz', 'meanRTmslo-hi40Hz', 'meanRTmslo-hi80Hz'))
myylim = c(min(Dtmp), max(Dtmp)) 
RmCI = tableCIs(Dtmp)
#RmCI
DRT = RmCI[1:3,1:4]

Dtmp = cbind(D$fp, Dtmp)
Dtmp = na.omit(Dtmp)
colnames(Dtmp) = c('fp','lo-hi20Hz','lo-hi40Hz','lo-hi80Hz')
bp2 = plotoverlay(dm = Dtmp, lbl = "RT", dvy = "RT (ms)", myylim)
grid.arrange(bp1, bp2, ncol = 2)

Ddpr$Condition = c("20 Hz", "40 Hz", "80 Hz")
DRT$Condition = c("20 Hz", "40 Hz", "80 Hz")
Ddpr$Condition = factor(Ddpr$Condition, levels=c("20 Hz", "40 Hz", "80 Hz"))
DRT$Condition = factor(DRT$Condition, levels=c("20 Hz", "40 Hz", "80 Hz"))
fig_dpr = ggplot(Ddpr, aes(x=Condition, y=Mean, fill=Condition)) + 
  geom_bar(stat="identity"
          ,color="black" # add black border to each bar
          ,position=position_dodge(), show.legend=FALSE) + # separare bars
  theme_bw() + # get rid of background
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank()) + #no grid lines
  geom_hline(yintercept = 0)  +
  geom_errorbar(position=position_dodge(.9), width=.25, 
                aes(ymin=LL, ymax=UL)) + 
  scale_fill_manual(
        values = c("#999999", "#CCCCCC", "#FFFFFF")) +
  labs(title = paste0("d´(N = ", nrow(D),")")) +
    theme(plot.title = element_text(hjust = 0.5)) +
  labs(x = "Low minus high load") +
  labs(y = "Sensitivity (d´)")
fig_RT = ggplot(DRT, aes(x=Condition, y=Mean, fill=Condition)) + 
  geom_bar(stat="identity"
          ,color="black" # add black border to each bar
          ,position=position_dodge(), show.legend=FALSE) + # separare bars
  theme_bw() + # get rid of background
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank()) + #no grid lines
  geom_hline(yintercept = 0)  +
  geom_errorbar(position=position_dodge(.9), width=.25, 
                aes(ymin=LL, ymax=UL)) + 
  scale_fill_manual(
        values = c("#999999", "#CCCCCC", "#FFFFFF")) +
  labs(title = paste0("RT (N = ", nrow(D),")")) +
    theme(plot.title = element_text(hjust = 0.5)) +
  labs(x = "Low minus high load") +
  labs(y = "Reaction time (ms)")
grid.arrange(fig_dpr, fig_RT, ncol = 2)

bp1 = arrangeGrob(fig_dpr, fig_RT, ncol = 2)
ggsave(file.path(dir_fig,"figure_performance.jpg"), bp1, dpi = 600, width = 7, height = 4)
ggsave(file.path(dir_fig,"figure_performance.pdf"), bp1, dpi = 600, width = 7, height = 4)
rm(Ddpr, DRT, fig_dpr, fig_RT, bp1, RmCI)

ANOVA

Analyze the factorial design (load x frequency) with repeated-measures ANOVA:

load: low and high
modulation frequency: 20, 40, 80 Hz

Dtmp1a = subset(D, select=c(fp, dprlo20Hz, dprlo40Hz, dprlo80Hz, 
                                dprhi20Hz, dprhi40Hz, dprhi80Hz))
Dtmp1b = subset(D, select=c(fp, meanRTmslo20Hz, meanRTmslo40Hz, meanRTmslo80Hz, 
                                meanRTmshi20Hz, meanRTmshi40Hz, meanRTmshi80Hz))
Dtmp1 = list(Dtmp1a, Dtmp1b)
Dvars = c("d´",'RT')
Fres = (matrix(ncol = 4, nrow = length(Dvars)))
Fres = as.data.frame(Fres)
colnames(Fres) = c('Variable','Overall_Interaction_p','BF01_Interaction', 'BF01_error%')
d prime
Fres = computeANOVA(1, Dtmp1[[1]], mylvls = c('low','high'))

## =====================================================================
##                         d´
## =====================================================================
## 
## 
## Descriptives
## ============
##  load frequency  N     Mean    CI_LL    CI_UL
##   low        20 25 4.390445 4.017362 4.763527
##   low        40 25 4.324651 4.046021 4.603281
##   low        80 25 4.225636 3.917314 4.533958
##  high        20 25 2.556847 2.253343 2.860351
##  high        40 25 2.628708 2.343701 2.913715
##  high        80 25 2.557600 2.295689 2.819511
## 
## 
## Anova
## =====
## 
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
## 
##                 Sum Sq num Df Error SS den Df F value    Pr(>F)    
## (Intercept)    1782.60      1   49.603     24 862.497 < 2.2e-16 ***
## load            112.56      1   10.073     24 268.191 1.581e-14 ***
## frequency         0.23      2    9.637     48   0.580    0.5638    
## load:frequency    0.20      2    8.697     48   0.542    0.5851    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Mauchly Tests for Sphericity
## 
##                Test statistic p-value
## frequency             0.99529 0.94717
## load:frequency        0.74031 0.03150
## 
## 
## Greenhouse-Geisser and Huynh-Feldt Corrections
##  for Departure from Sphericity
## 
##                 GG eps Pr(>F[GG])
## frequency      0.99531     0.5630
## load:frequency 0.79385     0.5458
## 
##                   HF eps Pr(>F[HF])
## frequency      1.0851207  0.5637804
## load:frequency 0.8408845  0.5555762
## 
## 
## BF without vs with interaction
## ==============================
## Bayes factor analysis
## --------------
## [1] load + frequency + fp : 6.279871 ±3.12%
## 
## Against denominator:
##   dv ~ load + frequency + load:frequency + fp 
## ---
## Bayes factor type: BFlinearModel, JZS
RT (ms)
Fres = computeANOVA(2, Dtmp1[[2]], mylvls = c('low','high'))

## =====================================================================
##                         RT
## =====================================================================
## 
## 
## Descriptives
## ============
##  load frequency  N     Mean    CI_LL    CI_UL
##   low        20 25 372.0967 349.2083 394.9851
##   low        40 25 374.7445 350.7671 398.7220
##   low        80 25 370.0598 352.2752 387.8444
##  high        20 25 486.6405 462.1193 511.1617
##  high        40 25 496.9289 469.9772 523.8806
##  high        80 25 494.4596 470.4060 518.5133
## 
## 
## Anova
## =====
## 
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
## 
##                  Sum Sq num Df Error SS den Df   F value    Pr(>F)    
## (Intercept)    28056925      1   358149     24 1880.1276 < 2.2e-16 ***
## load             543389      1    78899     24  165.2908  2.96e-12 ***
## frequency          1050      2    14363     48    1.7542    0.1840    
## load:frequency      668      2    16419     48    0.9771    0.3838    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Mauchly Tests for Sphericity
## 
##                Test statistic p-value
## frequency             0.87231 0.20782
## load:frequency        0.96035 0.62795
## 
## 
## Greenhouse-Geisser and Huynh-Feldt Corrections
##  for Departure from Sphericity
## 
##                 GG eps Pr(>F[GG])
## frequency      0.88676     0.1883
## load:frequency 0.96186     0.3812
## 
##                   HF eps Pr(>F[HF])
## frequency      0.9524285  0.1858081
## load:frequency 1.0439485  0.3837704
## 
## 
## BF without vs with interaction
## ==============================
## Bayes factor analysis
## --------------
## [1] load + frequency + fp : 6.734905 ±2.3%
## 
## Against denominator:
##   dv ~ load + frequency + load:frequency + fp 
## ---
## Bayes factor type: BFlinearModel, JZS

ANOVA summary

Analyze the factorial design (load x frequency) with repeated-measures ANOVA:

load: low and high
modulation frequency: 20, 40, 80 Hz

An interaction would suggest that load effects differed among the modulation frequencies.

Summary of interaction analyses:
* Overall_interaction has 2 dfs.
* BF01_interaction is the BF for a model without the overall interaction versus a model with the overall interaction (prior: Cauchy r scale = 0.5).
* BF01_error% estimates the error of the BF in percent. For example, if BF = 10 and the error is 10%, then the BF is estimated to vary between 9 and 11.

Fres[,2:4] = round(Fres[,2:4], digits = 3)
Fres[Fres < 0.001] = 0.001
Fres %>%
  kable(digits = 3) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), 
                full_width = F, position = 'left')
Variable Overall_Interaction_p BF01_Interaction BF01_error%
d´ 0.546 6.280 3.115
RT 0.381 6.735 2.303

visual P3

visual P3 = mean amplitude (between 300 and 400 ms) to targets minus nontargets
Note: Because there was no task during no load, there were no targets and nontargets to examine the P3.

table

Dtmp = subset(D, select=c('vP3lo20', 'vP3hi20', 'vP3lo-hi20',
                          'vP3lo40', 'vP3hi40', 'vP3lo-hi40',
                          'vP3lo80', 'vP3hi80', 'vP3lo-hi80'))
RmCI = tableCIs(Dtmp)
RmCI %>%
  kable(digits = 3) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), 
                full_width = F, position = 'left')
Variable Mean LL UL N
vP3lo20 9.696 7.371 12.021 25
vP3hi20 2.934 1.459 4.409 25
vP3lo-hi20 6.761 4.684 8.839 25
vP3lo40 9.685 7.914 11.457 25
vP3hi40 3.095 1.494 4.696 25
vP3lo-hi40 6.591 5.256 7.925 25
vP3lo80 9.957 7.971 11.943 25
vP3hi80 3.738 2.410 5.066 25
vP3lo-hi80 6.219 4.410 8.028 25

figure

Dtmp = subset(D, select=c('vP3lo-hi20', 'vP3lo-hi40', 'vP3lo-hi80'))
RmCI = tableCIs(Dtmp)
myylim = c(min(Dtmp), max(Dtmp)) 
Dtmp = cbind(D$fp, Dtmp)
Dtmp = na.omit(Dtmp)
colnames(Dtmp) = c('fp','lo-hi20Hz','lo-hi40Hz','lo-hi80Hz')
bp1 = plotoverlay(dm = Dtmp,lbl = "Visual P3", dvy = "Mean amplitude (µV)", myylim)
plot(bp1)

RmCI$Variable = c("20 Hz", "40 Hz", "80 Hz")
fig_P3 = ggplot(RmCI, aes(x=Variable, y=Mean, fill=Variable)) + 
  geom_bar(stat="identity"
          ,color="black" # add black border to each bar
          ,position=position_dodge(), show.legend=TRUE) + # separare bars
  theme_bw() + # get rid of background
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank()) + #no grid lines
  geom_hline(yintercept = 0)  +
  geom_errorbar(position=position_dodge(.9), width=.25, 
                aes(ymin=LL, ymax=UL)) + 
  scale_fill_manual(
        values = c("#999999", "#CCCCCC", "#FFFFFF"),
        name = "Frequency") +
  labs(title = paste0("Visual P3 (N = ", nrow(D),")")) +
    theme(plot.title = element_text(hjust = 0.5)) +
  labs(x = "Low minus high load") +
  labs(y = "Mean amplitude (µV)")
plot(fig_P3)

filename = "figure_visual_P3"
ggsave(file.path(dir_fig,paste0(filename,'.jpg')), fig_P3, dpi = 600, width = 7, height = 4)
ggsave(file.path(dir_fig,paste0(filename,'.pdf')), fig_P3, dpi = 600, width = 7, height = 4)
#rm(fig_P3, RmCI, filename)

ANOVA

Analyze the factorial design (load x frequency) with repeated-measures ANOVA:

load: low and high
modulation frequency: 20, 40, 80 Hz

Dtmp1a = subset(D, select=c(fp, vP3lo20, vP3lo40, vP3lo80, 
                                vP3hi20, vP3hi40, vP3hi80))
# rename because Hz is missing at the end
colnames(Dtmp1a) = c('fp', 'vP3lo20Hz', 'vP3lo40Hz', 'vP3lo80Hz', 
                           'vP3hi20Hz', 'vP3hi40Hz', 'vP3hi80Hz')
Dtmp1 = list(Dtmp1a)
Dvars = c("Visual P3")
Fres = (matrix(ncol = 4, nrow = length(Dvars)))
Fres = as.data.frame(Fres)
colnames(Fres) = c('Variable','Overall_Interaction_p','BF01_Interaction', 'BF01_error%')
Fres = computeANOVA(1, Dtmp1[[1]], mylvls = c('low','high'))

## =====================================================================
##                         Visual P3
## =====================================================================
## 
## 
## Descriptives
## ============
##  load frequency  N     Mean    CI_LL     CI_UL
##   low        20 25 9.695603 7.370569 12.020637
##   low        40 25 9.685139 7.913580 11.456698
##   low        80 25 9.957019 7.971284 11.942755
##  high        20 25 2.934121 1.458943  4.409298
##  high        40 25 3.094605 1.493526  4.695685
##  high        80 25 3.737903 2.410278  5.065529
## 
## 
## Anova
## =====
## 
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
## 
##                Sum Sq num Df Error SS den Df F value    Pr(>F)    
## (Intercept)    6371.5      1  1863.13     24 82.0747 3.262e-09 ***
## load           1596.0      1   422.25     24 90.7109 1.263e-09 ***
## frequency         8.3      2   152.01     48  1.3122    0.2787    
## load:frequency    1.9      2   237.40     48  0.1943    0.8240    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Mauchly Tests for Sphericity
## 
##                Test statistic  p-value
## frequency             0.68232 0.012327
## load:frequency        0.89686 0.285989
## 
## 
## Greenhouse-Geisser and Huynh-Feldt Corrections
##  for Departure from Sphericity
## 
##                 GG eps Pr(>F[GG])
## frequency      0.75891     0.2752
## load:frequency 0.90651     0.8031
## 
##                   HF eps Pr(>F[HF])
## frequency      0.7994202  0.2761494
## load:frequency 0.9763668  0.8190030
## 
## 
## BF without vs with interaction
## ==============================
## Bayes factor analysis
## --------------
## [1] load + frequency + fp : 7.785318 ±2.52%
## 
## Against denominator:
##   dv ~ load + frequency + load:frequency + fp 
## ---
## Bayes factor type: BFlinearModel, JZS

ANOVA summary

Analyze the factorial design (load x frequency) with repeated-measures ANOVA:

load: low and high
modulation frequency: 20, 40, 80 Hz

An interaction would suggest that load effects differed among the modulation frequencies.

Summary of interaction analyses:
* Overall_interaction has 2 dfs.
* BF01_interaction is the BF for a model without the overall interaction versus a model with the overall interaction (prior: Cauchy r scale = 0.5).
* BF01_error% estimates the error of the BF in percent. For example, if BF = 10 and the error is 10%, then the BF is estimated to vary between 9 and 11.

Fres[,2:4] = round(Fres[,2:4], digits = 3)
Fres[Fres < 0.001] = 0.001
Fres %>%
  kable(digits = 3) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), 
                full_width = F, position = 'left')
Variable Overall_Interaction_p BF01_Interaction BF01_error%
Visual P3 0.803 7.785 2.522

workload ratings (NASA-TLX)

workload ratings = mental demand and effort

table

Dnasa = data.frame()
tmpRmCI = data.frame(row.names = NULL)[1:18, ]
for (c in c('Mental', 'Effort')){
  Dtmp = getdiffs(D, iv = c)
  Dtmp = Dtmp[,-c(1:9)]
  myylim = c(min(Dtmp), max(Dtmp)) 
  RmCI = tableCIs(Dtmp)
  #write.table(RmCI, file = paste0("Table_descr_",c,"_",nowsample,"_subjects.txt"),
  #            row.names = F, dec = ".", sep = "\t")
  tmpRmCI = cbind(tmpRmCI, RmCI)
  
  Dtmp2 = RmCI
  Dtmp2$var = c
  Dnasa = rbind(Dnasa, Dtmp2)
  
  Dtmp = cbind(D$fp, Dtmp)
  Dtmp = na.omit(Dtmp)
  colnames(Dtmp)[1] = 'fp'
  #bp1 = plotoverlay(dm = Dtmp, lbl = paste0("NASA ", c), dvy = "Mean rating", myylim)
  #plot(bp1)
}
rownames(tmpRmCI) = NULL
tmpRmCI = tmpRmCI[,-c(5,6)]
tmpRmCI = subset(tmpRmCI, select=c(Variable, N, Mean, LL, UL, Mean.1, LL.1, UL.1))
tmpRmCI %>%
  kable(digits = 3,
        #caption = paste0("Workload ratings (N = ", nrow(D),")"),
        col.names = c('Condition','N','Mean','LL','UL','Mean','LL','UL')) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = F, position = 'left') %>%
  add_header_above(c(" " = 2, "Mental demand" = 3, "Effort" = 3))
Mental demand
Effort
Condition N Mean LL UL Mean LL UL
no20Hz 25 18.00 11.305 24.695 17.00 9.912 24.088
lo20Hz 25 28.44 21.885 34.995 29.50 22.728 36.272
hi20Hz 25 40.64 34.081 47.199 43.40 36.793 50.007
no-lo20Hz 25 -10.44 -16.938 -3.942 -12.50 -18.938 -6.062
no-hi20Hz 25 -22.64 -29.833 -15.447 -26.40 -33.942 -18.858
lo-hi20Hz 25 -12.20 -16.615 -7.785 -13.90 -20.450 -7.350
no40Hz 25 16.22 10.480 21.960 15.20 9.512 20.888
lo40Hz 25 26.62 20.774 32.466 27.94 20.214 35.666
hi40Hz 25 40.94 34.455 47.425 44.94 38.974 50.906
no-lo40Hz 25 -10.40 -15.603 -5.197 -12.74 -19.443 -6.037
no-hi40Hz 25 -24.72 -31.144 -18.296 -29.74 -35.224 -24.256
lo-hi40Hz 25 -14.32 -17.542 -11.098 -17.00 -22.421 -11.579
no80Hz 25 16.12 10.531 21.709 15.38 10.273 20.487
lo80Hz 25 29.28 23.230 35.330 32.26 24.136 40.384
hi80Hz 25 41.00 34.200 47.800 43.92 37.086 50.754
no-lo80Hz 25 -13.16 -18.154 -8.166 -16.88 -22.520 -11.240
no-hi80Hz 25 -24.88 -32.633 -17.127 -28.54 -35.421 -21.659
lo-hi80Hz 25 -11.72 -16.219 -7.221 -11.66 -17.232 -6.088
colnames(Dnasa)[1] = 'Condition'
colnames(Dnasa)[6] = 'Variable'
tmp = Dnasa$Condition %in% c('no20Hz', 'lo20Hz', 'hi20Hz',
                             'no40Hz', 'lo40Hz', 'hi40Hz',
                             'no80Hz', 'lo80Hz', 'hi80Hz')
Dnasa = Dnasa[tmp,]
Dnasa$Condition = factor(Dnasa$Condition, 
                         levels=c('no20Hz', 'lo20Hz', 'hi20Hz',
                                  'no40Hz', 'lo40Hz', 'hi40Hz',
                                  'no80Hz', 'lo80Hz', 'hi80Hz'))
Dnasa$Variable = factor(Dnasa$Variable, levels=c("Mental", "Effort"))
rm(c, bp1, bp2, tmp, Dtmp, Dtmp2, tmpRmCI, RmCI)

figure

colnames(Dnasa)[1] = 'Cond'
Dnasa$Freq = NA
Dnasa$Freq[grepl(pattern = "20", Dnasa$Cond)] = 20
Dnasa$Freq[grepl(pattern = "40", Dnasa$Cond)] = 40
Dnasa$Freq[grepl(pattern = "80", Dnasa$Cond)] = 80
Dnasa$Freq = factor(Dnasa$Freq)
Dnasa$Condition = NA
Dnasa$Condition[grepl(pattern = "no", Dnasa$Cond)] = 'no'
Dnasa$Condition[grepl(pattern = "lo", Dnasa$Cond)] = 'low'
Dnasa$Condition[grepl(pattern = "hi", Dnasa$Cond)] = 'high'
Dnasa$Condition = factor(Dnasa$Condition, 
                         levels=c('no', 'low', 'high'))
Dtmp = Dnasa[Dnasa$Variable=='Mental',]
fig_mental = ggplot(Dtmp, aes(x=Freq, y=Mean, fill=Condition)) + 
  geom_bar(stat="identity"
           ,color="black" # add black border to each bar
           ,position=position_dodge(), show.legend=TRUE) + # separate bars
  theme_bw() + # get rid of background
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank()) + #no grid lines
  geom_hline(yintercept = 0)  +
  geom_errorbar(position=position_dodge(.9), width=.25, 
                aes(ymin=LL, ymax=UL)) + 
  scale_fill_manual(
    values = c("#FFFFFF", "#CCCCCC", "#999999")) +
  labs(title = paste0("Mental demand (N = ", nrow(D),")")) +
  theme(plot.title = element_text(hjust = 0.5)) +
  labs(x = "Frequency (Hz)") +
  labs(y = "Mean rating") + 
  ylim(0, 50)
Dtmp = Dnasa[Dnasa$Variable=='Effort',]
fig_effort = ggplot(Dtmp, aes(x=Freq, y=Mean, fill=Condition)) + 
  geom_bar(stat="identity"
           ,color="black" # add black border to each bar
           ,position=position_dodge(), show.legend=TRUE) + # separate bars
  theme_bw() + # get rid of background
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank()) + #no grid lines
  geom_hline(yintercept = 0)  +
  geom_errorbar(position=position_dodge(.9), width=.25, 
                aes(ymin=LL, ymax=UL)) + 
  scale_fill_manual(
    values = c("#FFFFFF", "#CCCCCC", "#999999")) +
  labs(title = paste0("Effort (N = ", nrow(D),")")) +
  theme(plot.title = element_text(hjust = 0.5)) +
  labs(x = "Frequency (Hz)") +
  labs(y = "Mean rating") +
  ylim(0, 50)
grid.arrange(fig_mental, fig_effort, ncol = 2)

bp1 = arrangeGrob(fig_mental, fig_effort, ncol = 2)
filename = "figure_NASA_ratings"
ggsave(file.path(dir_fig,paste0(filename,'.jpg')), bp1, dpi = 600, width = 7, height = 4)
ggsave(file.path(dir_fig,paste0(filename,'.pdf')), bp1, dpi = 600, width = 7, height = 4)
rm(bp1, Dnasa, Dtmp, fig_mental, fig_effort, filename)

ANOVA

Analyze the factorial design (load x frequency) with repeated-measures ANOVA:

load: no, low, high
modulation frequency: 20, 40, 80 Hz

Dtmp1a = subset(D, select=c(fp, Mentalno20Hz, Mentalno40Hz, Mentalno80Hz, 
                                Mentallo20Hz, Mentallo40Hz, Mentallo80Hz,
                                Mentalhi20Hz, Mentalhi40Hz, Mentalhi80Hz))
Dtmp1b = subset(D, select=c(fp, Effortno20Hz, Effortno40Hz, Effortno80Hz, 
                                Effortlo20Hz, Effortlo40Hz, Effortlo80Hz,  
                                Efforthi20Hz, Efforthi40Hz, Efforthi80Hz))
Dtmp1 = list(Dtmp1a, Dtmp1b)
Dvars = c("Mental", "Effort")
Fres = (matrix(ncol = 4, nrow = length(Dvars)))
Fres = as.data.frame(Fres)
colnames(Fres) = c('Variable','Overall_Interaction_p','BF01_Interaction', 'BF01_error%')
Mental demand
Fres = computeANOVA(1, Dtmp1[[1]], mylvls = c('no','low','high'))

## =====================================================================
##                         Mental
## =====================================================================
## 
## 
## Descriptives
## ============
##  load frequency  N  Mean    CI_LL    CI_UL
##    no        20 25 18.00 11.30531 24.69469
##    no        40 25 16.22 10.47989 21.96011
##    no        80 25 16.12 10.53053 21.70947
##   low        20 25 28.44 21.88480 34.99520
##   low        40 25 26.62 20.77443 32.46557
##   low        80 25 29.28 23.23008 35.32992
##  high        20 25 40.64 34.08069 47.19931
##  high        40 25 40.94 34.45503 47.42497
##  high        80 25 41.00 34.19957 47.80043
## 
## 
## Anova
## =====
## 
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
## 
##                Sum Sq num Df Error SS den Df  F value    Pr(>F)    
## (Intercept)    183841      1    34889     24 126.4629 4.735e-11 ***
## load            21769      2     9907     48  52.7388 7.668e-13 ***
## frequency          51      2     1019     48   1.1920    0.3124    
## load:frequency    100      4     4059     96   0.5891    0.6713    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Mauchly Tests for Sphericity
## 
##                Test statistic p-value
## load                  0.33634  0.0000
## frequency             0.92532  0.4096
## load:frequency        0.43112  0.0268
## 
## 
## Greenhouse-Geisser and Huynh-Feldt Corrections
##  for Departure from Sphericity
## 
##                 GG eps Pr(>F[GG])    
## load           0.60109  1.379e-08 ***
## frequency      0.93051     0.3106    
## load:frequency 0.73639     0.6211    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                   HF eps   Pr(>F[HF])
## load           0.6152843 9.714385e-09
## frequency      1.0055900 3.124417e-01
## load:frequency 0.8506436 6.446848e-01
## 
## 
## BF without vs with interaction
## ==============================
## Bayes factor analysis
## --------------
## [1] load + frequency + fp : 24.42255 ±6.95%
## 
## Against denominator:
##   dv ~ load + frequency + load:frequency + fp 
## ---
## Bayes factor type: BFlinearModel, JZS
Effort
Fres = computeANOVA(2, Dtmp1[[2]], mylvls = c('no','low','high'))

## =====================================================================
##                         Effort
## =====================================================================
## 
## 
## Descriptives
## ============
##  load frequency  N  Mean     CI_LL    CI_UL
##    no        20 25 17.00  9.912279 24.08772
##    no        40 25 15.20  9.512408 20.88759
##    no        80 25 15.38 10.273404 20.48660
##   low        20 25 29.50 22.727550 36.27245
##   low        40 25 27.94 20.213619 35.66638
##   low        80 25 32.26 24.135708 40.38429
##  high        20 25 43.40 36.793375 50.00662
##  high        40 25 44.94 38.973905 50.90609
##  high        80 25 43.92 37.086064 50.75394
## 
## 
## Anova
## =====
## 
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
## 
##                Sum Sq num Df Error SS den Df  F value    Pr(>F)    
## (Intercept)    201811      1    37808     24 128.1064 4.152e-11 ***
## load            29878      2    11345     48  63.2087 3.561e-14 ***
## frequency          50      2     2697     48   0.4493    0.6407    
## load:frequency    269      4     5362     96   1.2023    0.3149    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Mauchly Tests for Sphericity
## 
##                Test statistic p-value
## load                  0.97116 0.71427
## frequency             0.93798 0.47889
## load:frequency        0.74933 0.69316
## 
## 
## Greenhouse-Geisser and Huynh-Feldt Corrections
##  for Departure from Sphericity
## 
##                 GG eps Pr(>F[GG])    
## load           0.97197  7.698e-14 ***
## frequency      0.94160     0.6292    
## load:frequency 0.88364     0.3155    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                  HF eps   Pr(>F[HF])
## load           1.056367 3.561477e-14
## frequency      1.019138 6.407296e-01
## load:frequency 1.055005 3.149139e-01
## 
## 
## BF without vs with interaction
## ==============================
## Bayes factor analysis
## --------------
## [1] load + frequency + fp : 13.56714 ±2.62%
## 
## Against denominator:
##   dv ~ load + frequency + load:frequency + fp 
## ---
## Bayes factor type: BFlinearModel, JZS

ANOVA summary

Analyze the factorial design (load x frequency) with repeated-measures ANOVA:

load: no, low, high
modulation frequency: 20, 40, 80 Hz

An interaction would suggest that load effects differed among the modulation frequencies.

Summary of interaction analyses:
* Overall_interaction has 4 dfs.
* BF01_interaction is the BF for a model without the overall interaction versus a model with the overall interaction (prior: Cauchy r scale = 0.5).
* BF01_error% estimates the error of the BF in percent. For example, if BF = 10 and the error is 10%, then the BF is estimated to vary between 9 and 11.

Fres[,2:4] = round(Fres[,2:4], digits = 3)
Fres[Fres < 0.001] = 0.001
Fres %>%
  kable(digits = 3) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), 
                full_width = F, position = 'left')
Variable Overall_Interaction_p BF01_Interaction BF01_error%
Mental 0.621 24.423 6.953
Effort 0.315 13.567 2.625

ASSR means and prereg BFs

Means and preregistered Bayes factors.

Ampl SmN

Amplitude signal minus noise (SmN)

plot

iv = "AmpSmN"
Dtmp = getdiffs(D, iv = iv)
myylim = c(min(Dtmp), max(Dtmp)) 
Dtmp = Dtmp[,c(10:27)]
Dtmp = cbind(D$fp, Dtmp)
colnames(Dtmp)[1] = 'fp'
Dtmp1 = Dtmp[,c(1,2:7)]
colnames(Dtmp1) = c('fp', 'no', 'low', 'high', 'no-low', 'no-high', 'low-high')
bp1a = plotoverlay(dm = Dtmp1,lbl = paste0("20 Hz"), dvy = "Mean SmN amplitude (µV)", myylim)
Dtmp1 = Dtmp[,c(1,8:13)]
colnames(Dtmp1) = c('fp', 'no', 'low', 'high', 'no-low', 'no-high', 'low-high')
bp1b = plotoverlay(dm = Dtmp1,lbl = paste0("40 Hz"), dvy = "", myylim)
Dtmp1 = Dtmp[,c(1,14:19)]
colnames(Dtmp1) = c('fp', 'no', 'low', 'high', 'no-low', 'no-high', 'low-high')
bp1c = plotoverlay(dm = Dtmp1,lbl = paste0("80 Hz"), dvy = "", myylim)
grid.arrange(bp1a, bp1b, bp1c, ncol = 3)

bp1 = arrangeGrob(bp1a, bp1b, bp1c, ncol = 3)
ggsave(file.path(dir_fig,"figure_amp.jpg"), bp1, dpi = 600, width = 7, height = 4)
ggsave(file.path(dir_fig,"figure_amp.pdf"), bp1, dpi = 600, width = 7, height = 4)
Dtmp_amp = Dtmp # save for later

table

iv = "AmpSmN"
Dtmp = getdiffs(D, iv = iv)
Dtmp = tableCIs(Dtmp[,-c(1:9)])
Dtmp %>%
  kable(digits = 3) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), 
                full_width = F, position = 'left')
Variable Mean LL UL N
no20Hz 0.360 0.263 0.457 25
lo20Hz 0.333 0.225 0.441 25
hi20Hz 0.348 0.236 0.460 25
no-lo20Hz 0.027 -0.028 0.082 25
no-hi20Hz 0.012 -0.056 0.080 25
lo-hi20Hz -0.015 -0.063 0.034 25
no40Hz 0.455 0.340 0.570 25
lo40Hz 0.420 0.316 0.524 25
hi40Hz 0.411 0.303 0.518 25
no-lo40Hz 0.035 0.009 0.061 25
no-hi40Hz 0.045 0.010 0.080 25
lo-hi40Hz 0.010 -0.027 0.046 25
no80Hz 0.089 0.047 0.131 25
lo80Hz 0.096 0.056 0.137 25
hi80Hz 0.088 0.043 0.133 25
no-lo80Hz -0.007 -0.031 0.016 25
no-hi80Hz 0.001 -0.023 0.026 25
lo-hi80Hz 0.008 -0.018 0.035 25

preregistered BFs

The BFs are computed from Bayesian one-sample t tests of difference scores (with Aladins R script).
https://doi.org/10.17045/sthlmuni.4981154.v3

The BF01 uses uniform H1 models with different lower limits (LL) and upper limits (UL).
BF01 is the evidence for the null hypothesis relative to the alternative hypothesis.
If BF01 > 3, this is evidence for the null. If BF01 < 1/3, this is evidence against the null.

Dtmp = contrastBFs(D, iv = iv)
Dtmp %>%
  kable(digits = 3,
        col.names = c('Variable','Mean','[-1, +1]','[0, +1]','[0, +0.24]','[0, +0.37]')) %>%
  add_header_above(c(" " = 2, "BF01" = 4)) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), 
                full_width = F, position = 'left')
BF01
Variable Mean [-1, +1] [0, +1] [0, +0.24] [0, +0.37]
no-lo20Hz 0.027 17.658 10.459 2.520 3.882
no-hi20Hz 0.012 22.409 17.290 4.177 6.432
lo-hi20Hz -0.015 27.508 50.083 12.254 18.829
no-lo40Hz 0.035 1.874 0.942 0.226 0.348
no-hi40Hz 0.045 1.948 0.981 0.235 0.363
lo-hi40Hz 0.010 38.429 26.962 6.537 10.061
no-lo80Hz -0.007 56.473 99.945 24.944 38.197
no-hi80Hz 0.001 65.824 58.949 14.471 22.223
lo-hi80Hz 0.008 49.496 32.773 7.963 12.251

ITC SmN

Intertrial coherence signal minus noise (SmN)

plot

iv = "ItcSmN"
Dtmp = getdiffs(D, iv = iv)
myylim = c(min(Dtmp), max(Dtmp)) 
Dtmp = Dtmp[,c(10:27)]
Dtmp = cbind(D$fp, Dtmp)
colnames(Dtmp)[1] = 'fp'
Dtmp1 = Dtmp[,c(1,2:7)]
colnames(Dtmp1) = c('fp', 'no', 'low', 'high', 'no-low', 'no-high', 'low-high')
bp1a = plotoverlay(dm = Dtmp1,lbl = paste0("20 Hz"), dvy = "Mean SmN ITC", myylim)
Dtmp1 = Dtmp[,c(1,8:13)]
colnames(Dtmp1) = c('fp', 'no', 'low', 'high', 'no-low', 'no-high', 'low-high')
bp1b = plotoverlay(dm = Dtmp1,lbl = paste0("40 Hz"), dvy = "", myylim)
Dtmp1 = Dtmp[,c(1,14:19)]
colnames(Dtmp1) = c('fp', 'no', 'low', 'high', 'no-low', 'no-high', 'low-high')
bp1c = plotoverlay(dm = Dtmp1,lbl = paste0("80 Hz"), dvy = "", myylim)
grid.arrange(bp1a, bp1b, bp1c, ncol = 3)

bp1 = arrangeGrob(bp1a, bp1b, bp1c, ncol = 3)
ggsave(file.path(dir_fig,"figure_itc.jpg"), bp1, dpi = 600, width = 7, height = 4)
ggsave(file.path(dir_fig,"figure_itc.pdf"), bp1, dpi = 600, width = 7, height = 4)
Dtmp_itc = Dtmp # save for later

table

iv = "ItcSmN"
Dtmp = getdiffs(D, iv = iv)
Dtmp = tableCIs(Dtmp[,-c(1:9)])
Dtmp %>%
  kable(digits = 3) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), 
                full_width = F, position = 'left')
Variable Mean LL UL N
no20Hz 0.172 0.122 0.222 25
lo20Hz 0.184 0.137 0.231 25
hi20Hz 0.181 0.125 0.238 25
no-lo20Hz -0.012 -0.030 0.006 25
no-hi20Hz -0.009 -0.039 0.021 25
lo-hi20Hz 0.003 -0.024 0.029 25
no40Hz 0.352 0.279 0.425 25
lo40Hz 0.326 0.257 0.395 25
hi40Hz 0.302 0.227 0.377 25
no-lo40Hz 0.026 -0.006 0.058 25
no-hi40Hz 0.050 0.022 0.079 25
lo-hi40Hz 0.024 -0.011 0.059 25
no80Hz 0.116 0.068 0.165 25
lo80Hz 0.104 0.054 0.154 25
hi80Hz 0.094 0.050 0.138 25
no-lo80Hz 0.012 -0.013 0.037 25
no-hi80Hz 0.023 -0.005 0.050 25
lo-hi80Hz 0.010 -0.019 0.039 25

preregistered BFs

The BFs are computed from Bayesian one-sample t tests of difference scores (with Aladins R script).
https://doi.org/10.17045/sthlmuni.4981154.v3

The BF01 uses uniform H1 models with different lower limits (LL) and upper limits (UL).
BF01 is the evidence for the null hypothesis relative to the alternative hypothesis.
If BF01 > 3, this is evidence for the null. If BF01 < 1/3, this is evidence against the null.

Dtmp = contrastBFs(D, iv = iv)
Dtmp %>%
  kable(digits = 3,
        col.names = c('Variable','Mean','[-1, +1]','[0, +1]','[0, +0.24]','[0, +0.37]')) %>%
  add_header_above(c(" " = 2, "BF01" = 4)) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), 
                full_width = F, position = 'left')
BF01
Variable Mean [-1, +1] [0, +1] [0, +0.24] [0, +0.37]
no-lo20Hz -0.012 37.180 171.635 44.156 67.259
no-hi20Hz -0.009 44.121 78.747 19.487 29.884
lo-hi20Hz 0.002 60.065 50.804 12.432 19.102
no-lo40Hz 0.026 12.555 6.602 1.588 2.448
no-hi40Hz 0.050 0.241 0.121 0.029 0.045
lo-hi40Hz 0.024 16.840 9.138 2.201 3.391
no-lo80Hz 0.012 38.449 22.652 5.483 8.440
no-hi80Hz 0.022 14.620 7.687 1.850 2.851
lo-hi80Hz 0.010 43.010 27.771 6.735 10.365

Floor effects?

Does the combined effect of task and load (i.e., no load minus high load) for SmN depend on the initial SmN value for no load?
A floor effect would be supported if large difference scores of no minus high load (on Y axis) can be observed only for large initial values of no load (on X axis).
In the plots, the regression lines are derived from local polynomial regression fitting, see geom_smooth(method=‘loess’).

Amp SmN

Results suggest no obvious relationship.
Note: 40-Hz shows an increase, but this is driven by two data points.

# amp SmN
Dtmp = Dtmp_amp
tymin = min(Dtmp$'no-hi20Hz',Dtmp$'no-hi40Hz',Dtmp$'no-hi80Hz')
tymax = max(Dtmp$'no-hi20Hz',Dtmp$'no-hi40Hz',Dtmp$'no-hi80Hz')
bp1 = ggplot(Dtmp, aes(x=`no20Hz`, y=`no-hi20Hz`)) +
  theme_bw() + # get rid of background
  geom_point(shape=1, size = 4) +    # Use hollow circles
  labs(title = paste0("20 Hz (N = ",nrow(D),")")) +
  labs(x = "No load SmN") +
  labs(y = "No minus high SmN") +
  geom_hline(yintercept = 0) +
  ylim(tymin,tymax) +
  geom_smooth(method = "loess", formula = y ~ x)
bp2 = ggplot(Dtmp, aes(x=`no40Hz`, y=`no-hi40Hz`)) +
  theme_bw() + # get rid of background
  geom_point(shape=1, size = 4) +    # Use hollow circles
  labs(title = paste0("40 Hz (N = ",nrow(D),")")) +
  labs(x = "No load SmN") +
  labs(y = "No minus high SmN") +
  geom_hline(yintercept = 0) +
  ylim(tymin,tymax) +
  geom_smooth(method = "loess", formula = y ~ x)
bp3 = ggplot(Dtmp, aes(x=`no80Hz`, y=`no-hi80Hz`)) +
  theme_bw() + # get rid of background
  geom_point(shape=1, size = 4) +    # Use hollow circles
  labs(title = paste0("80 Hz (N = ",nrow(D),")")) +
  labs(x = "No load SmN") +
  labs(y = "No minus high SmN") +
  geom_hline(yintercept = 0) +
  ylim(tymin,tymax) +
  geom_smooth(method = "loess", formula = y ~ x)
grid.arrange(bp1, bp2, bp3, ncol = 3)

bp3 = arrangeGrob(bp1, bp2, ncol = 2)
ggsave(file.path(dir_fig,"figure_scatter_amp.jpg"), bp3, dpi = 600, width = 13, height = 4)
ggsave(file.path(dir_fig,"figure_scatter_amp.pdf"), bp3, dpi = 600, width = 13, height = 4)
rm(tymax, tymin, bp1, bp2, bp3)

ITC SmN

Results suggest no obvious relationship.
Note: 80-Hz shows an increase, but this is driven by few data points.

# itc SmN
Dtmp = Dtmp_itc
tymin = min(Dtmp$'no-hi20Hz',Dtmp$'no-hi40Hz',Dtmp$'no-hi80Hz')
tymax = max(Dtmp$'no-hi20Hz',Dtmp$'no-hi40Hz',Dtmp$'no-hi80Hz')
bp1 = ggplot(Dtmp, aes(x=`no20Hz`, y=`no-hi20Hz`)) +
  theme_bw() + # get rid of background
  geom_point(shape=1, size = 4) +    # Use hollow circles
  labs(title = paste0("20 Hz (N = ",nrow(D),")")) +
  labs(x = "No load SmN") +
  labs(y = "No minus high SmN") +
  geom_hline(yintercept = 0) +
  ylim(tymin,tymax) +
  geom_smooth(method = "loess", formula = y ~ x)
bp2 = ggplot(Dtmp, aes(x=`no40Hz`, y=`no-hi40Hz`)) +
  theme_bw() + # get rid of background
  geom_point(shape=1, size = 4) +    # Use hollow circles
  labs(title = paste0("40 Hz (N = ",nrow(D),")")) +
  labs(x = "No load SmN") +
  labs(y = "No minus high SmN") +
  geom_hline(yintercept = 0) +
  ylim(tymin,tymax) +
  geom_smooth(method = "loess", formula = y ~ x)
bp3 = ggplot(Dtmp, aes(x=`no80Hz`, y=`no-hi80Hz`)) +
  theme_bw() + # get rid of background
  geom_point(shape=1, size = 4) +    # Use hollow circles
  labs(title = paste0("80 Hz (N = ",nrow(D),")")) +
  labs(x = "No load SmN") +
  labs(y = "No minus high SmN") +
  geom_hline(yintercept = 0) +
  ylim(tymin,tymax) +
  geom_smooth(method = "loess", formula = y ~ x)
grid.arrange(bp1, bp2, bp3, ncol = 3)

bp3 = arrangeGrob(bp1, bp2, ncol = 2)
ggsave(file.path(dir_fig,"figure_scatter_itc.jpg"), bp3, dpi = 600, width = 13, height = 4)
ggsave(file.path(dir_fig,"figure_scatter_itc.pdf"), bp3, dpi = 600, width = 13, height = 4)
rm(tymax, tymin, bp1, bp2, bp3, Dtmp)

Exploratory BFs

In the preregistered BF analyses, we defined upper limits of the priors on the basis of the mean SmN values for no load in a previous study with 40-Hz ASSRs (Study 2 in Szychowska & Wiens, 2020). Theoretically, the SmN during no load defines an upper limit of how strong the task effect and load effect can be.

Because the visual task in the present study was identical to that in the previous study, we anticipated that we would observe similar mean SmN values during no load for 40-Hz. So, we used our previous data to predict the upper limits for 40 Hz. For simplicity, we used these upper limits also for 20 and 80 Hz. Thus, the original goal was to define realistic upper limits in the BF analyses and to preregister these limits to facilitate confirmatory research (hypothesis testing).

However, the present results showed that the actual SmN during no and low load differed from the preregistered upper limits. To conduct BF analyses that matched our original goal of using reasonable upper limits, we defined the upper limits on the basis of the present data, separately for each measure and frequency.

No minus low/high
One analysis considered no load as the baseline, and the H1 model was defined as a uniform distribution between 0 (as lower limit, LL) and the present mean of no load (as upper limit, UL).

No minus low = task effect
No minus high = combined effect of task and load

Low minus high
The other analysis considered low load as the baseline, and the H1 model was defined as a uniform distribution between 0 (as lower limit, LL) and the present mean of low load (as upper limit, UL).

Low minus high = load effect

In the table, LL and UL refer to the 95% CIs of the means (and not to the LL and UL of the H1 models).

no minus low/high

# create matrix with measure, freq, mean no load, 
# diff noMlo, LL_CI, UL_CI, BF01 noMlo, 
# diff noMhi, LL_CI, UL_CI, BF01 noMhi
DBF01 = numeric()
Dtmp_freq = c(20,40,80)

# amp ----
Dtmp_H1 = subset(Dtmp_amp, 
                 select=c(no20Hz, no40Hz, no80Hz))
Dtmp_noMlo = subset(Dtmp_amp, 
                    select=c('no-lo20Hz', 'no-lo40Hz', 'no-lo80Hz'))
Dtmp_noMhi = subset(Dtmp_amp, 
                    select=c('no-hi20Hz', 'no-hi40Hz', 'no-hi80Hz'))
LL = 0
for (i in 1:3) {
  # define UL as mean of no load
  UL = mean(Dtmp_H1[,i])
  
  # noMlo
  tmp = Dtmp_noMlo[,i]
  # noMhi
  tmp2 = Dtmp_noMhi[,i]
  
  DBF01 = rbind(DBF01, c(Dtmp_freq[i], UL, 
                         mean(tmp),  ci(tmp),  BFfitH1(LL = 0, UL, tmp),
                         mean(tmp2), ci(tmp2), BFfitH1(LL = 0, UL, tmp2)))
}
rm(Dtmp_H1, Dtmp_noMlo, Dtmp_noMhi)

# itc ----
Dtmp_H1 = subset(Dtmp_itc, 
                 select=c(no20Hz, no40Hz, no80Hz))
Dtmp_noMlo = subset(Dtmp_itc, 
                    select=c('no-lo20Hz', 'no-lo40Hz', 'no-lo80Hz'))
Dtmp_noMhi = subset(Dtmp_itc, 
                    select=c('no-hi20Hz', 'no-hi40Hz', 'no-hi80Hz'))
LL = 0
for (i in 1:3) {
  # define UL as mean of no load
  UL = mean(Dtmp_H1[,i])
  
  # noMlo
  tmp = Dtmp_noMlo[,i]
  # noMhi
  tmp2 = Dtmp_noMhi[,i]
  
  DBF01 = rbind(DBF01, c(Dtmp_freq[i], UL, 
                         mean(tmp),  ci(tmp),  BFfitH1(LL = 0, UL, tmp),
                         mean(tmp2), ci(tmp2), BFfitH1(LL = 0, UL, tmp2)))
}
DBF01 = as.data.frame(DBF01)
colnames(DBF01)=c('Freq','M-no',
                  'M-noMlo','LL-noMlo','UL-noMlo','BF01-noMlo',
                  'M-noMhi','LL-noMhi','UL-noMhi','BF01-noMhi')
DBF01$Measure = c('Amp','Amp','Amp','ITC','ITC','ITC')
DBF01 = DBF01[,c(ncol(DBF01),1:ncol(DBF01)-1)]
DBF01 %>%
  kable(digits = 3,
        col.names = c('Measure','Frequency','Mean',
                      'Mean','LL', 'UL', 'BF01',
                      'Mean','LL', 'UL', 'BF01')) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = F, position = 'left') %>%
  add_header_above(c(" " = 2, "no" = 1, "no minus low" = 4, "no minus high" = 4))
no
no minus low
no minus high
Measure Frequency Mean Mean LL UL BF01 Mean LL UL BF01
Amp 20 0.360 0.027 -0.028 0.082 3.772 0.012 -0.056 0.080 6.250
Amp 40 0.455 0.035 0.009 0.061 0.428 0.045 0.010 0.080 0.446
Amp 80 0.089 -0.007 -0.031 0.016 9.336 0.001 -0.023 0.026 5.398
ITC 20 0.172 -0.012 -0.030 0.006 31.859 -0.009 -0.039 0.021 14.010
ITC 40 0.352 0.026 -0.006 0.058 2.329 0.050 0.022 0.079 0.043
ITC 80 0.116 0.012 -0.013 0.037 2.658 0.023 -0.005 0.050 0.896
rm(DBF01, Dtmp_freq, Dtmp_H1, Dtmp_noMlo, Dtmp_noMhi)

low minus high

# create matrix with measure, freq, mean low load, 
# diff loMhi, LL_CI, UL_CI, BF01 looMhi
DBF01 = numeric()
Dtmp_freq = c(20,40,80)

# amp ----
Dtmp_H1 = subset(Dtmp_amp, 
                 select=c(lo20Hz, lo40Hz, lo80Hz))
Dtmp_loMhi = subset(Dtmp_amp, 
                    select=c('lo-hi20Hz', 'lo-hi40Hz', 'lo-hi80Hz'))
LL = 0
for (i in 1:3) {
  # define UL as mean of no load
  UL = mean(Dtmp_H1[,i])
  
  # loMhi
  tmp = Dtmp_loMhi[,i]
  
  DBF01 = rbind(DBF01, c(Dtmp_freq[i], UL, 
                         mean(tmp),  ci(tmp),  BFfitH1(LL = 0, UL, tmp)))
}
rm(Dtmp_H1, Dtmp_loMhi)

# itc ----
Dtmp_H1 = subset(Dtmp_itc, 
                 select=c(lo20Hz, lo40Hz, lo80Hz))
Dtmp_loMhi = subset(Dtmp_itc, 
                    select=c('lo-hi20Hz', 'lo-hi40Hz', 'lo-hi80Hz'))
LL = 0
for (i in 1:3) {
  # define UL as mean of no load
  UL = mean(Dtmp_H1[,i])
  
  # loMhi
  tmp = Dtmp_loMhi[,i]
  
  DBF01 = rbind(DBF01, c(Dtmp_freq[i], UL, 
                         mean(tmp),  ci(tmp),  BFfitH1(LL = 0, UL, tmp)))
}
DBF01 = as.data.frame(DBF01)
colnames(DBF01)=c('Freq','M-no',
                  'M-loMhi','LL-loMhi','UL-loMhi','BF01-loMhi')
DBF01$Measure = c('Amp','Amp','Amp','ITC','ITC','ITC')
DBF01 = DBF01[,c(ncol(DBF01),1:ncol(DBF01)-1)]
DBF01 %>%
  kable(digits = 3,
        col.names = c('Measure','Frequency','Mean',
                      'Mean','LL', 'UL', 'BF01')) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = F, position = 'left') %>%
  add_header_above(c(" " = 2, "low" = 1, "low minus high" = 4))
low
low minus high
Measure Frequency Mean Mean LL UL BF01
Amp 20 0.333 -0.015 -0.063 0.034 16.961
Amp 40 0.420 0.010 -0.027 0.046 11.404
Amp 80 0.096 0.008 -0.018 0.035 3.202
ITC 20 0.184 0.003 -0.024 0.029 9.531
ITC 40 0.326 0.024 -0.011 0.059 2.987
ITC 80 0.104 0.010 -0.019 0.039 2.920
rm(DBF01, Dtmp_freq, Dtmp_H1, Dtmp_loMhi)

40-Hz: combine with previous data

Our previous study (Szychowska & Wiens, 2020) studied effects of task and load on 40-Hz ASSRs.

Exp. 1: Experiment 1 (N = 43) in Szychowska & Wiens (2020). This study included low and high load.

Exp. 2: Experiment 2 (N = 45) in Szychowska & Wiens (2020). This study included no, low, high, and very high load. Because results showed that performance was impaired more strongly during high load than very high load, high load was used below (and very high load was excluded).

Below, we combined the previous data (Exp. 1 and 2) with the present data to obtain an estimate on the basis of all of our data on 40-Hz ASSRs.

Dtmpa = subset(Dtmp_amp, 
               select=c(fp, no40Hz, lo40Hz, hi40Hz))
colnames(Dtmpa) = c('fp', 'AmpSmNno', 'AmpSmNlo', 'AmpSmNhi')
Dtmpi = subset(Dtmp_itc, 
               select=c(no40Hz, lo40Hz, hi40Hz))
colnames(Dtmpi) = c('ItcSmNno', 'ItcSmNlo', 'ItcSmNhi')
Dtmp = cbind(Dtmpa, Dtmpi)
Dtmp$study = 'present'

DS1=read.csv(file.path(dir_res,"data_export_SmN_ASSR_study1_all_subjects.tsv"), sep="\t",header=TRUE)
DS2=read.csv(file.path(dir_res,"data_export_SmN_ASSR_study2.tsv"), sep="\t",header=TRUE)

no minus low/high

No load is used as the upper limit of the H1 model.
H1 model:
LL = 0
UL = no load
Note that Exp. 1 is excluded because no load was not administered.

No minus low = task effect
No minus high = combined effect of task and load

In the table, LL and UL refer to the 95% CIs of the means (and not to the LL and UL of the H1 models).

# compare noMlo and noMhi

Dall = merge(Dtmp, DS2, all = T)

DBF01 = numeric()

# create matrix with measure, study, n, mean no load, 
# diff noMlo, LL 95%CI, UL 95%CI, BF01 noMlo,
# diff noMhi, LL 95%CI, UL 95%CI, BF01 noMhi
LL = 0

# amp
# study 2
UL = mean(DS2$AmpSmNno)
tmp = DS2$AmpSmNno-DS2$AmpSmNlo
tmp2 = DS2$AmpSmNno-DS2$AmpSmNhi
DBF01 = rbind(DBF01, c('Amp','Exp. 2', nrow(DS2), UL, 
                       mean(tmp), ci(tmp), BFfitH1(LL = 0, UL, tmp),
                       mean(tmp2), ci(tmp2), BFfitH1(LL = 0, UL, tmp2)))
# present study
UL = mean(Dtmp$AmpSmNno)
tmp = Dtmp$AmpSmNno-Dtmp$AmpSmNlo
tmp2 = Dtmp$AmpSmNno-Dtmp$AmpSmNhi
DBF01 = rbind(DBF01, c('Amp','Present', nrow(Dtmp), UL, 
                       mean(tmp), ci(tmp), BFfitH1(LL = 0, UL, tmp),
                       mean(tmp2), ci(tmp2), BFfitH1(LL = 0, UL, tmp2)))
# combined
UL = mean(Dall$AmpSmNno)
tmp = Dall$AmpSmNno-Dall$AmpSmNlo
tmp2 = Dall$AmpSmNno-Dall$AmpSmNhi
DBF01 = rbind(DBF01, c('Amp','Combined', nrow(Dall), UL, 
                       mean(tmp), ci(tmp), BFfitH1(LL = 0, UL, tmp),
                       mean(tmp2), ci(tmp2), BFfitH1(LL = 0, UL, tmp2)))

# itc
# study 2
UL = mean(DS2$ItcSmNno)
tmp = DS2$ItcSmNno-DS2$ItcSmNlo
tmp2 = DS2$ItcSmNno-DS2$ItcSmNhi
DBF01 = rbind(DBF01, c('ITC','Exp. 2', nrow(DS2), UL, 
                       mean(tmp), ci(tmp), BFfitH1(LL = 0, UL, tmp),
                       mean(tmp2), ci(tmp2), BFfitH1(LL = 0, UL, tmp2)))
# present study
UL = mean(Dtmp$ItcSmNno)
tmp = Dtmp$ItcSmNno-Dtmp$ItcSmNlo
tmp2 = Dtmp$ItcSmNno-Dtmp$ItcSmNhi
DBF01 = rbind(DBF01, c('ITC','Present', nrow(Dtmp), UL, 
                       mean(tmp), ci(tmp), BFfitH1(LL = 0, UL, tmp),
                       mean(tmp2), ci(tmp2), BFfitH1(LL = 0, UL, tmp2)))
# combined
UL = mean(Dall$ItcSmNno)
tmp = Dall$ItcSmNno-Dall$ItcSmNlo
tmp2 = Dall$ItcSmNno-Dall$ItcSmNhi
DBF01 = rbind(DBF01, c('ITC','Combined', nrow(Dall), UL, 
                       mean(tmp), ci(tmp), BFfitH1(LL = 0, UL, tmp),
                       mean(tmp2), ci(tmp2), BFfitH1(LL = 0, UL, tmp2)))
DBF01 = as.data.frame(DBF01)
colnames(DBF01)=c('Measure','Study','N','Mno',
                  'M-noMlo','LL-noMlo','UL-noMlo','BF01-noMlo',
                  'M-noMhi','LL-noMhi','UL-noMhi','BF01-noMhi')
# convert character to numeric
i = 3:12
DBF01[, i] = apply(DBF01[,i],2, function(x) as.numeric(as.character(x)))
DBF01 %>%
  kable(digits = 3,
        col.names = c('Measure','Study','N','Mean',
                      'Mean','LL', 'UL', 'BF01',
                      'Mean','LL', 'UL', 'BF01')) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = F, position = 'left') %>%
  add_header_above(c(" " = 3, "no" = 1, "no minus low" = 4, "no minus high" = 4)) %>%
  row_spec(c(3,6), bold = T)
no
no minus low
no minus high
Measure Study N Mean Mean LL UL BF01 Mean LL UL BF01
Amp Exp. 2 45 0.243 -0.001 -0.012 0.010 39.213 0.000 -0.013 0.012 32.045
Amp Present 25 0.455 0.035 0.009 0.061 0.428 0.045 0.010 0.080 0.446
Amp Combined 70 0.319 0.012 0.000 0.024 3.039 0.016 0.000 0.031 2.133
ITC Exp. 2 45 0.368 -0.004 -0.021 0.012 51.978 0.005 -0.015 0.024 19.749
ITC Present 25 0.352 0.026 -0.006 0.058 2.329 0.050 0.022 0.079 0.043
ITC Combined 70 0.362 0.006 -0.009 0.022 16.239 0.021 0.004 0.037 0.820

low minus high

Low load is used as the upper limit of the H1 model.
H1 model:
LL = 0
UL = low load

Low minus high = load effect

In the table, LL and UL refer to the 95% CIs of the means (and not to the LL and UL of the H1 models).

# compare low with high load
Dall = merge(Dtmp, DS1, all = T)
Dall = merge(Dall, DS2, all = T)

DBF01 = numeric()

# create matrix with measure, study, n, mean low load, 
# diff loMhi, LL 95%CI, UL 95%CI, BF01 loMhi
LL = 0

# amp
# study 1
UL = mean(DS1$AmpSmNlo)
tmp = DS1$AmpSmNlo-DS1$AmpSmNhi
DBF01 = rbind(DBF01, c('Amp','Exp. 1', nrow(DS1), UL, mean(tmp), ci(tmp), BFfitH1(LL = 0, UL, tmp)))
# study 2
UL = mean(DS2$AmpSmNlo)
tmp = DS2$AmpSmNlo-DS2$AmpSmNhi
DBF01 = rbind(DBF01, c('Amp','Exp. 2', nrow(DS2), UL, mean(tmp), ci(tmp), BFfitH1(LL = 0, UL, tmp)))
# present study
UL = mean(Dtmp$AmpSmNlo)
tmp = Dtmp$AmpSmNlo-Dtmp$AmpSmNhi
DBF01 = rbind(DBF01, c('Amp','Present', nrow(Dtmp), UL, mean(tmp), ci(tmp), BFfitH1(LL = 0, UL, tmp)))
# combined
UL = mean(Dall$AmpSmNlo)
tmp = Dall$AmpSmNlo-Dall$AmpSmNhi
DBF01 = rbind(DBF01, c('Amp','Combined', nrow(Dall), UL, mean(tmp), ci(tmp), BFfitH1(LL = 0, UL, tmp)))

# itc
# study 1
UL = mean(DS1$ItcSmNlo)
tmp = DS1$ItcSmNlo-DS1$ItcSmNhi
DBF01 = rbind(DBF01, c('ITC','Exp. 1', nrow(DS1), UL, mean(tmp), ci(tmp), BFfitH1(LL = 0, UL, tmp)))
# study 2
UL = mean(DS2$ItcSmNlo)
tmp = DS2$ItcSmNlo-DS2$ItcSmNhi
DBF01 = rbind(DBF01, c('ITC','Exp. 2', nrow(DS2), UL, mean(tmp), ci(tmp), BFfitH1(LL = 0, UL, tmp)))
# present study
UL = mean(Dtmp$ItcSmNlo)
tmp = Dtmp$ItcSmNlo-Dtmp$ItcSmNhi
DBF01 = rbind(DBF01, c('ITC','Present', nrow(Dtmp), UL, mean(tmp), ci(tmp), BFfitH1(LL = 0, UL, tmp)))
# combined
UL = mean(Dall$ItcSmNlo)
tmp = Dall$ItcSmNlo-Dall$ItcSmNhi
DBF01 = rbind(DBF01, c('ITC','Combined', nrow(Dall), UL, mean(tmp), ci(tmp), BFfitH1(LL = 0, UL, tmp)))
DBF01 = as.data.frame(DBF01)
colnames(DBF01)=c('Measure','Study','N','Mean','Mean','LL-95%CI', 'UL-95%CI', 'BF01')
# convert character to numeric
i = 3:8
DBF01[, i] = apply(DBF01[,i],2, function(x) as.numeric(as.character(x)))
DBF01 %>%
  kable(digits = 3,
        col.names = c('Measure','Study','N','Mean','Mean','LL', 'UL', 'BF01')) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = F, position = 'left') %>%
  add_header_above(c(" " = 3, "low" = 1, "low minus high" = 4)) %>%
  row_spec(c(4,8), bold = T)
low
low minus high
Measure Study N Mean Mean LL UL BF01
Amp Exp. 1 43 0.191 0.006 -0.004 0.017 7.760
Amp Exp. 2 45 0.244 0.000 -0.012 0.013 28.960
Amp Present 25 0.420 0.010 -0.027 0.046 11.404
Amp Combined 113 0.262 0.005 -0.005 0.015 16.130
ITC Exp. 1 43 0.309 0.019 0.000 0.038 1.924
ITC Exp. 2 45 0.373 0.009 -0.009 0.027 11.840
ITC Present 25 0.326 0.024 -0.011 0.059 2.987
ITC Combined 113 0.338 0.016 0.004 0.028 0.879