Project website

Data and other resources: Open Science Framework

  
## here() starts at /Volumes/psychology$/Researchers/reteig/AB_tDCS-sEBR
## ── Attaching packages ──────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.2.1 ──
## ✔ ggplot2 3.1.0       ✔ purrr   0.2.5  
## ✔ tibble  2.1.1       ✔ dplyr   0.8.0.1
## ✔ tidyr   0.8.3       ✔ stringr 1.3.1  
## ✔ readr   1.1.1       ✔ forcats 0.3.0
## Warning: package 'tibble' was built under R version 3.5.2
## Warning: package 'tidyr' was built under R version 3.5.2
## Warning: package 'dplyr' was built under R version 3.5.2
## ── Conflicts ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## 
## Attaching package: 'psych'
## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha
## Loading required package: coda
## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## The following object is masked from 'package:tidyr':
## 
##     expand
## ************
## Welcome to BayesFactor 0.9.12-4.2. If you have questions, please contact Richard Morey (richarddmorey@gmail.com).
## 
## Type BFManual() to open the manual.
## ************
## Loading required package: nlme
## 
## Attaching package: 'nlme'
## The following object is masked from 'package:dplyr':
## 
##     collapse
## This is mgcv 1.8-24. For overview type 'help("mgcv-package")'.
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## Loading required package: hypergeo
## R version 3.5.1 (2018-07-02)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: OS X El Capitan 10.11.6
## 
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/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] hypergeo_1.2-13        lmtest_0.9-36          zoo_1.8-4             
##  [4] sandwich_2.5-0         mgcv_1.8-24            nlme_3.1-137          
##  [7] cocor_1.1-3            BayesFactor_0.9.12-4.2 Matrix_1.2-14         
## [10] coda_0.19-2            broom_0.5.1            psych_1.8.10          
## [13] forcats_0.3.0          stringr_1.3.1          dplyr_0.8.0.1         
## [16] purrr_0.2.5            readr_1.1.1            tidyr_0.8.3           
## [19] tibble_2.1.1           ggplot2_3.1.0          tidyverse_1.2.1       
## [22] here_0.1               rmdformats_0.3.3       knitr_1.21            
## 
## loaded via a namespace (and not attached):
##  [1] httr_1.3.1         jsonlite_1.6       elliptic_1.3-9    
##  [4] modelr_0.1.2       gtools_3.8.1       shiny_1.2.0       
##  [7] assertthat_0.2.0   highr_0.7          cellranger_1.1.0  
## [10] yaml_2.2.0         pillar_1.3.1       backports_1.1.3   
## [13] lattice_0.20-35    glue_1.3.0         digest_0.6.18     
## [16] promises_1.0.1     rvest_0.3.2        colorspace_1.3-2  
## [19] htmltools_0.3.6    httpuv_1.4.5       plyr_1.8.4        
## [22] pkgconfig_2.0.2    haven_1.1.2        questionr_0.6.3   
## [25] bookdown_0.9       xtable_1.8-4       mvtnorm_1.0-8     
## [28] scales_1.0.0       later_0.7.5        MatrixModels_0.4-1
## [31] generics_0.0.2     withr_2.1.2        pbapply_1.3-4     
## [34] lazyeval_0.2.1     cli_1.0.1          mnormt_1.5-5      
## [37] magrittr_1.5       crayon_1.3.4       readxl_1.1.0      
## [40] mime_0.6           evaluate_0.12      MASS_7.3-50       
## [43] xml2_1.2.0         foreign_0.8-70     tools_3.5.1       
## [46] hms_0.4.2          munsell_0.5.0      compiler_3.5.1    
## [49] contfrac_1.1-12    rlang_0.3.4        grid_3.5.1        
## [52] rstudioapi_0.8     miniUI_0.1.1.1     rmarkdown_1.11    
## [55] gtable_0.2.0       deSolve_1.21       R6_2.3.0          
## [58] lubridate_1.7.4    rprojroot_1.3-2    stringi_1.2.4     
## [61] parallel_3.5.1     Rcpp_1.0.0         tidyselect_0.2.5  
## [64] xfun_0.4

Prepare data

Demographics

## Using ',' as decimal and '.' as grouping mark. Use read_delim() for more control.
participant demographics and questionnaire data
subject session polarity time gender age sleep_hours sleep_quality alcohol_glasses smoker caffeine_today caffeine_usual contacts_glasses EEG_experience
S01 1 anodal 2015-02-20 09:00:00 female 20 7.5 4 0 no 0 2.0 yes no
S01 2 cathodal 2015-02-27 09:00:00 female 20 7.5 5 0 no 0 2.0 yes no
S02 1 cathodal 2015-02-22 14:00:00 female 19 8.0 2 2 no 1 1.0 no no
S02 2 anodal 2015-03-01 14:00:00 female 19 9.0 2 2 no 2 3.0 no no
S03 1 cathodal 2015-02-24 09:00:00 female 24 6.0 3 4 yes 0 4.5 no no
S03 2 anodal 2015-03-03 09:00:00 female 24 6.5 4 0 yes 2 3.0 no no
  1. subject: Participant ID, e.g. S01, S12
  2. session: Whether data is from the first (1) or 2nd (2) session
  3. polarity: Whether participant received anodal or cathodal tDCS
  4. time: Date-time the session was scheduled to start. Original format is DAY DD.MM.YYYY HH:MM (e.g. “Fri 10.06.2016 14:00”); parsed here as YYYY-MM-DD HH:MM:SS.
  5. gender
  6. age: in years
  7. sleep_hours: amount of sleep gotten the night before, in hours
  8. sleep_quality: quality of sleep, rated 1 (bad), 2 (subpar), 3 (OK), 4 (good), or 5 (very good)
  9. alcohol_glasses: number of alcoholic drinks consumed day before the session (or day of the session)
  10. smoker: whether participant considers themselves a smoker, yes or no
  11. caffeine_today: number of cups of coffee/black tea consumed on day of the session
  12. caffeine_usual: number of cups of coffee/black tea participant would typically consume
  13. contacts_glasses: whether participants wears contacts, glasses or neither (no). Unfortunately some (n = 6) also simply answered yes, without specifying. This data was collected because participants with contacts might have a higher spontaneous eye blink rate. However, even when counting all the yes responses as contacts-wearers, this would probably leave too few too properly analyze (n = 11).
  14. EEG_experience: whether participant previously participated in an EEG study, yes or no. This data was collected because participants with EEG experience might remember the instruction that blinks are bad for the signal, and would thus have a lower spontaneous eye blink rate. Only a handful (n=6) participants answered yes here though.

Reliability

sEBR

Plot

Plot session 1 data against session 2:

Reliability of spontaneous Eye Blink Rate

Reliability of spontaneous Eye Blink Rate

lower sEBR with EEG experience
EEG_experience n median_sEBR
no 34 13.3
yes 6 5.2

S27 could be classified as an outlier, because their sEBR in session 2 is the highest overall (57), and it is quite a bit higher in session 2 than in session 1 (difference of 12). This is not a data quality issue, as the blink component had little noise, and all the blinks seem to have been marked correctly.1

Statistics

Compute intra-class correlations (ICCs) for sEBR in session 1 and 2:

## Warning in if (class(x1) == as.character("try-error")) {: the condition has
## length > 1 and only the first element will be used
Intra-class correlations for sEBR
type ICC F df1 df2 p lower bound upper bound
Single_raters_absolute ICC1 0.845 11.937 39 40 0 0.728 0.915
Single_random_raters ICC2 0.845 11.683 39 39 0 0.726 0.915
Single_fixed_raters ICC3 0.842 11.683 39 39 0 0.721 0.913
Average_raters_absolute ICC1k 0.916 11.937 39 40 0 0.842 0.956
Average_random_raters ICC2k 0.916 11.683 39 39 0 0.841 0.956
Average_fixed_raters ICC3k 0.914 11.683 39 39 0 0.838 0.955

This computes all the ICCs from Shrout and Fleiss (1979).

Since this is a reliability/test-retest analysis, we want the single-rating, absolute agreement, two-way mixed effects model (or two-way random effects; same formula). This is ICC(2,1) in the conventions from Shrout and Fleiss (1979). See Koo & Li (2016) for more information.

The ICC is based on an ANOVA model, so we can also examine whether sEBR differs significantly between the sessions (and subjects, but that’s surely the case). There is no significant difference between sessions in sEBR: F(1, 39) = 0.149, p = 0.701. This is good, as at first we were worried participants would blink less in the second session, because they learned during the 1st session that blinks created artefacts in the EEG signal during task performance.

no difference in sEBR between sessions
session median_sEBR
1 11.7
2 12.6

The intra-class correlation is 0.85, with a 95% confidence interval of 0.73–0.92, indicating the level of reliability is “good” (Koo & Li, 2016).

Kruis et al. (2016) report Cronbach’s alpha values over three sessions for sEBR of 0.85 (long-term meditators) and 0.79 (meditation-naive participants). This is equivalent to an average-measures two-way ICC for consistency. In our data, this ICC is 0.91, 95%CI [0.84, 0.95]

AB magnitude

Statistics

Compute intra-class correlations for AB magnitude in session 1 and 2:

## Warning in if (class(x1) == as.character("try-error")) {: the condition has
## length > 1 and only the first element will be used
Intra-class correlations for AB magnitude
type ICC F df1 df2 p lower bound upper bound
Single_raters_absolute ICC1 0.574 3.694 39 40 0 0.325 0.749
Single_random_raters ICC2 0.598 5.128 39 39 0 0.246 0.790
Single_fixed_raters ICC3 0.674 5.128 39 39 0 0.461 0.813
Average_raters_absolute ICC1k 0.729 3.694 39 40 0 0.491 0.856
Average_random_raters ICC2k 0.748 5.128 39 39 0 0.394 0.882
Average_fixed_raters ICC3k 0.805 5.128 39 39 0 0.631 0.897

There is a significant difference in AB magnitude between the sessions: F(1, 39) = 16.53, p = 0. This can also be seen in the scatter plot: most participants have a smaller AB magnitude in the 2nd session.

The difference in AB magnitude over sessions seems to be driven by the short lag, so the AB genuinely becomes smaller in the second session. Usually this is not the case (only T2 performance improves over time, but no interaction with lag, Dale, Dux & Arnell (2013)).

The intra-class correlation is 0.6, indicating “moderate” reliability, with a 95% confidence interval of 0.25 (poor) – 0.79 (good).

For comparison with previous studies, let’s also look at the regular Pearson correlation:

## 
##  Pearson's product-moment correlation
## 
## data:  .$session_1 and .$session_2
## t = 5.7259, df = 38, p-value = 1.354e-06
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.4683432 0.8185404
## sample estimates:
##      cor 
## 0.680563

This is actually in the higher range of what is typically reported: (around .4–.6, Dale, Dux & Arnell (2013)).

Correlation sEBR ~ AB magnitude

Correlations between the two measures are bounded by their reliability, but since both were good (or at least acceptable), this sets a reasonable upper bound even with our relatively limited sample size.

Linear

Colzato et al. (2009) report a negative association between sEBR and AB magnitude: participants with higher sEBR have a smaller attentional blink. However, this relationship was not replicated by Georgopoulou & Slagter (2013).

Pearson

sEBR vs. AB magnitude correlation in either session
session estimate statistic p.value parameter conf.low conf.high method alternative
1 0.08 0.48 0.64 38 -0.24 0.38 Pearson’s product-moment correlation two.sided
2 0.00 0.02 0.99 38 -0.31 0.31 Pearson’s product-moment correlation two.sided

So here there’s virtually no association (even slightly positive in session 1).

Bayesian

One-sided

All prior mass folded to negative effect sizes only:

## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character
## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector
Bayes factor (in favor of H0) with one-sided prior
session rowname BF_null
1 Alt., r=0.333 -1<rho<0 3.90
2 Alt., r=0.333 -1<rho<0 2.87

Now the evidence for the null in the first session is a bit stronger (because there the direction of the effect is slightly positive).

Replication BF

Finally, we can specify a replication Bayes Factor [(Wagenmakers, Verhagen & Ly, 2016)]((https://doi.org/10.3758/s13428-015-0593-0) for the result in Colzato et al. (2008), who found a negative correlation: r(20) = -.530.2

Replication Bayes Factor (in favor of H0)
session bf0R
1 15.80
2 10.43

There is strong evidence in favor of the null hypothesis that there is no effect, vs. the hypothesis that the effect is as in Colzato et al. (2008).

“Inverted-u-shaped”

However, the relationship between dopamine levels (indexed by sEBR) and cognitive performance is usually not characterized as linear, but inverted u-shaped (Cools & d’Esposito, 2011).

This can be examined by testing for a sign-flip in the relationship between sEBR and AB magnitude: if it is inverted u-shaped, low-sEBR values (first half of the \(\cap\)) should correlate positively with ABmagnitude, whereas high-sEBR values (second half of the \(\cap\)) should correlate negatively. The two-lines test (Simonsohn, 2019) implements this by estimating a “break point” (the apex of the \(\cap\)) from the data, and then performing two regressions.

For the first session data, the overall shape is an inverted-U, but neither slope is significant (although the positive slope is trending).

The 2nd-session data looks different: there are more break-points in the cubic spline fit, and both slopes are (a tiny bit) positive. Still both-slopes are non-significant.

Correlation sEBR ~ tDCS effect on AB magnitude

Correlate AB magnitude change scores

The effects of tDCS could depend on an individual’s striatal dopamine levels. sEBR can be used as a proxy for the latter. For the former, we’ll use the “change from baseline” scores in the anodal and cathodal sessions, by subtracting baseline T2|T1 accuracy from T2|T1 accuracy during the tDCS or post blocks:

Pearson

sEBR vs. effect of tDCS correlations
polarity change estimate statistic p.value parameter conf.low conf.high method alternative
anodal tDCS - baseline 0.21 1.31 0.20 38 -0.11 0.49 Pearson’s product-moment correlation two.sided
anodal post - baseline 0.24 1.53 0.13 38 -0.08 0.51 Pearson’s product-moment correlation two.sided
cathodal tDCS - baseline -0.16 -1.02 0.31 38 -0.45 0.16 Pearson’s product-moment correlation two.sided
cathodal post - baseline 0.04 0.24 0.81 38 -0.28 0.35 Pearson’s product-moment correlation two.sided

Comparing the correlations in different blocks

A complementary way of looking at this is not to compute change scores, but to compute the correlation between sEBR and AB magnitude in each block.

sEBR vs. AB magnitude correlation per block
polarity block estimate statistic p.value parameter conf.low conf.high method alternative
anodal post 0.25 1.58 0.12 38 -0.07 0.52 Pearson’s product-moment correlation two.sided
anodal pre 0.04 0.23 0.82 38 -0.28 0.34 Pearson’s product-moment correlation two.sided
anodal tDCS 0.15 0.94 0.35 38 -0.17 0.44 Pearson’s product-moment correlation two.sided
cathodal post 0.09 0.53 0.60 38 -0.23 0.39 Pearson’s product-moment correlation two.sided
cathodal pre 0.05 0.29 0.78 38 -0.27 0.35 Pearson’s product-moment correlation two.sided
cathodal tDCS -0.05 -0.33 0.75 38 -0.36 0.26 Pearson’s product-moment correlation two.sided

The correlations can then also be directly compared to see if they differ significantly (e.g., whether the correlation in the “pre” block in the anodal/cathodal session is significantly different than the correlation in the “tDCS” or “post blocks”).

There are lots of tests for comparing correlations implemented in the cocor package; see its accompanying paper: https://doi.org/10.1371/journal.pone.0121945.

First we need the reliability for the AB magnitude scores, i.e. the correlation between the AB magnitude scores in the different blocks:

Reliability of AB magnitude in pre vs. subsequent blocks
polarity block estimate statistic p.value parameter conf.low conf.high method alternative
anodal tDCS 0.84 9.40 0 38 0.71 0.91 Pearson’s product-moment correlation two.sided
cathodal tDCS 0.82 8.78 0 38 0.68 0.90 Pearson’s product-moment correlation two.sided
anodal post 0.63 4.99 0 38 0.40 0.79 Pearson’s product-moment correlation two.sided
cathodal post 0.71 6.15 0 38 0.51 0.83 Pearson’s product-moment correlation two.sided

Now make a single table out of the 1) sEBR ~ AB magnitude correlations in each session/block 2) the pre-tDCS/pre-post reliabilities, where:

  • The r_baseline column has the correlation between AB magnitude and sEBR in the baseline block, for each polarity
  • The r_subsequent column houses the same correlation for the tDCS and post blocks (indicated by the block column)
  • The r_reliability column houses the correlation between AB magnitude in the pre block and one of the subsequent blocks (tDCS or post, indicated by the block column)
## Joining, by = c("polarity", "block")
Data frame for comparing correlations
polarity parameter r_baseline block r_subsequent r_reliability
anodal 38 0.04 tDCS 0.15 0.84
cathodal 38 0.05 tDCS -0.05 0.82
anodal 38 0.04 post 0.25 0.63
cathodal 38 0.05 post 0.09 0.71

Now we have all the information to perform a significance test of the 4 correlation pairs:

To fully characterize the correlation between sEBR and AB magnitude change scores, in addition to the three columns in the data frame, we’d also have to look at the difference in varability in AB magnitude scores between the blocks. See also Gardner & Neufeld (1987) and Griffin, Murray, & Gonzalez (1999).

Standard deviations for AB magnitude per block, as well as the ratio between the pre- and subsequent blocks (M)
polarity post pre tDCS M_tDCS M_post
anodal 0.1365154 0.1395991 0.1479759 1.0600060 0.9779101
cathodal 0.1195923 0.1362417 0.1281921 0.9409165 0.8777948

  1. We also tried rerunning the correlations without S27, but that did not qualitatively change the results (to do that, simply add filter(subject != "S27) %>% to the 2nd line of any of the code chunks).↩

  2. (I suppose it should be r(18), as it’s based on 20 data points)↩