Behavioral Data

Data Cleaning

---

# Load data
load(file = "./data/Single_Trial_Data.rda")


# Exclude missing responses, RT outliers, and trials with ERP artifacts
single_trial_data_clean <- single_trial_data %>%
  dplyr::filter(
      response_type != "miss" &
      rt_invalid    == FALSE &
      !is.na(MFN_0_100_FCz)
  ) # (53093 of 53760 trials left)


# Create numeric accuracy variable (1 = correct, 0 = incorrect)
single_trial_data_clean <- single_trial_data_clean %>%
  dplyr::mutate(accuracy_numeric = ifelse(response_type == "correct", 1, 0))


# Create column with single-trial PES (RT post-error - RT pre-error for all CCEC sequences)
single_trial_data_clean$pes <- NA
for (i in 3:(nrow(single_trial_data_clean) - 1)) {
  if (single_trial_data_clean[i,       ]$response_type == "incorrect" &
      single_trial_data_clean[(i + 1), ]$response_type == "correct" &
      single_trial_data_clean[(i - 1), ]$response_type == "correct" &
      single_trial_data_clean[(i - 2), ]$response_type == "correct" &
      # Exclude sequences in which a trial was excluded in between
      single_trial_data_clean[(i + 1), ]$trial - single_trial_data_clean[(i - 1), ]$trial == 2) {
    single_trial_data_clean[i, ]$pes <- (single_trial_data_clean[(i + 1), ]$rt) - (single_trial_data_clean[(i - 1), ]$rt)
  }
}


# For each last first and trial in a block, PES cannot be determined; set these values to NA
single_trial_data_clean[single_trial_data_clean$trial == 1   |
  single_trial_data_clean$trial == 80  | single_trial_data_clean$trial == 81 |
  single_trial_data_clean$trial == 160 | single_trial_data_clean$trial == 161 |
  single_trial_data_clean$trial == 240 | single_trial_data_clean$trial == 241 |
  single_trial_data_clean$trial == 320 | single_trial_data_clean$trial == 321 |
  single_trial_data_clean$trial == 400 | single_trial_data_clean$trial == 401 |
  single_trial_data_clean$trial == 480, "pes"] <- NA


# Make categorical variables factors
single_trial_data_clean$participant_id <- factor(single_trial_data_clean$participant_id)
single_trial_data_clean$group          <- factor(single_trial_data_clean$group)
single_trial_data_clean$stimulation    <- factor(single_trial_data_clean$stimulation)
single_trial_data_clean$response_type  <- factor(single_trial_data_clean$response_type)

Trials with a response time below 100 ms or above 800 ms or with a missing response were excluded from all analyses. We further discarded trials containing artifacts in the EEG, that is, a voltage change exceeding 50 μV between sample points or 200 μV within an epoch. Please see section “ERP Data” for the percentage of excluded trials.

Descriptive Statistics

---

This table corresponds to Table 2 in the manuscript.

# Calculate descriptive statistics for RT
descriptive_statistics_rt <- summarySEwithinO(
  data          = single_trial_data_clean,
  measurevar    = "rt",
  withinvars    = c("response_type", "stimulation"),
  betweenvars   = "group",
  idvar         = "participant_id",
  conf.interval = .95
) %>% dplyr::rename(M = rt, M_norm = rt_norm) %>%
  dplyr::arrange(response_type) %>%
  dplyr::select(-response_type)


# Calculate descriptive statistics for PES
descriptive_statistics_pes <- summarySEwithinO(
  data          = single_trial_data_clean[!is.na(single_trial_data_clean$pes), ],
  measurevar    = "pes",
  withinvars    = "stimulation",
  betweenvars   = "group",
  idvar         = "participant_id",
  conf.interval = .95
) %>% dplyr::rename(M = pes, M_norm = pes_norm)


# Calculate descriptive statistics for accuracy
descriptive_statistics_accuracy <- summarySEwithinO(
  data          = single_trial_data_clean,
  measurevar    = "accuracy_numeric",
  withinvars    = "stimulation",
  betweenvars   = "group",
  idvar         = "participant_id",
  conf.interval = .95
) %>% dplyr::rename(M = accuracy_numeric, M_norm = accuracy_numeric_norm) %>%
  # Multiply numeric values by 100 to obtain values in percent
  dplyr::mutate_if(is.numeric, list(~ . * 100))


# Combine measures
descriptive_statistics <- as.data.frame(rbind(descriptive_statistics_rt,
  descriptive_statistics_pes, descriptive_statistics_accuracy))


# Label measures
descriptive_statistics$Measure <- c(rep("RT correct (ms)", 4), rep("RT error (ms)", 4), 
                                    rep("PES (ms)", 4), rep("Accuracy (%)", 4))


# Create column with format "M [CI]" for RT and PES (round to 0 decimals)
descriptive_statistics$M_CI <- paste0(round(descriptive_statistics$M, 0), " [",
  round(descriptive_statistics$M - descriptive_statistics$ci, 0), ", ",
  round(descriptive_statistics$M + descriptive_statistics$ci, 0), "]")


# Create column with format "M [CI]" for accuracy (round to 2 decimals)
descriptive_statistics[c(13:16), ]$M_CI <- paste0(
  format(round(descriptive_statistics[c(13:16), ]$M,  2), nsmall = 2), " [",
  format(round(descriptive_statistics[c(13:16), ]$M -
               descriptive_statistics[c(13:16), ]$ci, 2), nsmall = 2), ", ",
  format(round(descriptive_statistics[c(13:16), ]$M +
               descriptive_statistics[c(13:16), ]$ci, 2), nsmall = 2), "]")


# Select relevant columns
descriptive_statistics <- descriptive_statistics[, c ("group", "stimulation", "Measure", "M_CI")]



# Split and re-merge table to display both groups next to each other
descriptive_statistics <- split(descriptive_statistics, list(descriptive_statistics$stimulation, 
                                                             descriptive_statistics$group))
descriptive_statistics_display <- cbind(
  descriptive_statistics$sham.OCD[, c("Measure", "M_CI")],
  descriptive_statistics$sham.HC[, "M_CI"],
  descriptive_statistics$verum.OCD[, "M_CI"],
  descriptive_statistics$verum.HC[, "M_CI"]
)


# Display table
my_table_template(descriptive_statistics_display,
  caption = "Behavioral Measures in Patients With OCD and HC Participants",
  col_names = c(" ", rep("M [95% CI]", 4)),
  header_above_config = c(" ", rep(c("Patients with OCD", "HC participants"), 2)),
  footnote_config = c(general = "CIs are adjusted for within-participant designs (Morey, 2008). 
  RT = response time; PES = post-error slowing.")) %>%
  add_header_above(c(" " = 1, "Sham tDCS" = 2, "Cathodal tDCS" = 2))

Behavioral Measures in Patients With OCD and HC Participants

	Sham tDCS		Cathodal tDCS
	Patients with OCD	HC participants	Patients with OCD	HC participants
	M [95% CI]	M [95% CI]	M [95% CI]	M [95% CI]
RT correct (ms)	416 [415, 418]	401 [400, 403]	411 [410, 412]	403 [401, 404]
RT error (ms)	346 [341, 352]	341 [335, 346]	345 [340, 350]	335 [330, 340]
PES (ms)	40 [28, 51]	29 [20, 39]	39 [29, 50]	32 [21, 42]
Accuracy (%)	94.68 [94.14, 95.21]	93.24 [92.64, 93.84]	94.14 [93.57, 94.70]	93.51 [92.92, 94.10]
Note:
CIs are adjusted for within-participant designs (Morey, 2008). RT = response time; PES = post-error slowing.

(G)LMM Analyses

---

We analyzed behavioral performance (response time, post-error slowing, and response accuracy) using (generalized) linear mixed models ([G]LMMs) on single-trial data. We tested whether group differences and tDCS effects were present. Group (healthy controls, OCD) and tDCS condition (cathodal, sham) were included as fixed effects in all models. All categorical fixed effects were effect-coded (contrast coefficients −0.5 and 0.5).

We determined the random-effects structure for each model based on the procedure proposed by Bates et al. (2015), starting with the maximal random-effects structure justified by the design, with by-participant random intercepts and random slopes for all fixed factors and (where applicable) their interactions. If required for model convergence, correlation parameters of the random terms were set to zero. Random effects preventing model convergence or explaining zero variance as determined by principal component analysis were removed to avoid overparameterization.

# Define contrasts (sliding difference contrasts = effect coding for factors with 2 levels)
contrasts(single_trial_data_clean$stimulation)   <- contr.sdif(2)
contrasts(single_trial_data_clean$group)         <- contr.sdif(2)
contrasts(single_trial_data_clean$response_type) <- contr.sdif(2)


# Prepare labels for (G)LMM tables
labels <- c(
  "(Intercept)"                               = "Intercept",
  "stimulation2-1"                            = "Stimulation",
  "group2-1"                                  = "Group",
  "response_type2-1"                          = "Response type",
  "stimulation2-1:group2-1"                   = "Stimulation x Group",
  "stimulation2-1:response_type2-1"           = "Stimulation x Response type",
  "group2-1:response_type2-1"                 = "Group x Response type",
  "stimulation2-1:group2-1:response_type2-1"  = "Stimulation x Group x Response type"
)

Response Time

Determine data transformation

Response time was log-transformed prior to analysis to meet the assumption of normally distributed residuals. The appropriate transformation was determined using the Box–Cox procedure (Box & Cox, 1964).

# Arrange plots
par(mfrow = c(1, 3))


# Determine transformation by estimating optimal lambda using Box–Cox procedure
bc_rt <- boxcox(rt ~ 1, data = single_trial_data_clean)
optlambda_rt <- bc_rt$x[which.max(bc_rt$y)]


# Density plot for log-transformed RT values
plot(density(single_trial_data_clean$rt_log), main = "Log RT: Density  Plot")


# Q-q plot for log-transformed RT values
qqnorm(single_trial_data_clean$rt_log, main = "Log RT: Q-Q Plot", pch = 1)

# Reset plot layout
par(mfrow = c(1, 1))

The optimal lambda is -0.38, suggesting that log transformation is appropriate.

LMM

We analyzed response time data using a LMM with response type (correct, error), group, and tDCS condition as predictors.

This table corresponds to Table S1 in the supplemental material.

# Run model
LMM_rt <- lmer(rt_log ~ stimulation * group * response_type +
  (1 + stimulation * response_type | participant_id),
data = single_trial_data_clean,
REML = TRUE,
control = lmerControl(optimizer = "bobyqa")
)


# Check model output and PCA of random-effects variance-covariance estimates
# summary(LMM_rt)        # Model does converge
# isSingular(LMM_rt)     # No singular fit
# summary(rePCA(LMM_rt)) # All terms explain variance


# Display results (fixed effects)
tab_model(LMM_rt,
  dv.labels = "log(RT)", pred.labels = labels, show.stat = TRUE, show.icc = FALSE, show.r2 = FALSE, 
  show.re.var = FALSE, show.ngroups = FALSE, minus.sign = "−", string.pred = "Fixed effects", 
  string.est = "b", string.stat = "t", string.ci = "95 % CI", p.val = "satterthwaite", 
  wrap.labels = 80, title = "Results of the LMM Predicting RT as a Function of Stimulation Condition 
  (Cathodal - Sham), Group (OCD - Healthy Controls), and Response Type (Incorrect - Correct)"
)


# Display random effects
print("Random effects:")
print(VarCorr(LMM_rt), digits = 1, comp = "Std.Dev.")

Results of the LMM Predicting RT as a Function of Stimulation Condition (Cathodal - Sham), Group (OCD - Healthy Controls), and Response Type (Incorrect - Correct)

	log(RT)
Fixed effects	b	95 % CI	t	p
Intercept	5.91	5.89 – 5.93	568.05	<0.001
Stimulation	-0.00	-0.02 – 0.01	-0.53	0.599
Group	0.02	-0.02 – 0.06	0.99	0.328
Response type	-0.17	-0.19 – -0.15	-16.46	<0.001
Stimulation x Group	0.00	-0.03 – 0.03	0.00	1.000
Stimulation x Response type	0.00	-0.02 – 0.02	0.14	0.890
Group x Response type	-0.02	-0.06 – 0.02	-0.80	0.426
Stimulation x Group x Response type	0.03	-0.01 – 0.07	1.68	0.099
Observations	53093

[1] "Random effects:"
 Groups         Name                            Std.Dev. Corr             
 participant_id (Intercept)                     0.08                      
                stimulation2-1                  0.05      0.04            
                response_type2-1                0.07      0.03  0.17      
                stimulation2-1:response_type2-1 0.06     -0.12  0.04 -0.18
 Residual                                       0.17                      

Response time analysis revealed a significant main effect of response type, indicating that responses were faster for errors than for correct responses. Response time did not differ significantly between groups or tDCS conditions and there was no significant interaction between any of the factors.

Post-Error Slowing

Determine data transformation

No data transformation was required to meet the assumption of normally distributed residuals. This was determined using the Box–Cox procedure (Box & Cox, 1964).

# Arrange plots
par(mfrow = c(1, 3))


# Determine transformation by estimating optimal lambda using Box–Cox procedure
bc_pes <- boxcox(pes + 1000 ~ 1, data = single_trial_data_clean[!is.na(single_trial_data_clean$pes), ])
optlambda_pes <- bc_pes$x[which.max(bc_pes$y)]


# Density plot for PES values
plot(density(single_trial_data_clean[!is.na(single_trial_data_clean$pes), ]$pes), main = "PES: Density  Plot")


# Q-q plot for PES values
qqnorm(single_trial_data_clean[!is.na(single_trial_data_clean$pes), ]$pes, main = "PES: Q-Q Plot", pch = 1)

# Reset plot layout
par(mfrow = c(1, 1))

The optimal lambda is 0.75, suggesting that no transformation is needed.

LMM

For the analysis of PES, we fitted a LMM on single-trial values that were calculated as the response time difference between correct responses that directly preceded and followed an error. This PES quantification results in a measure that is not confounded by fluctuations in motivation or response caution over time (Dutilh et al., 2012). In this analysis, we considered only error trials that were preceded by at least two correct trials and followed by at least one correct trial. Model estimates of the LMM on PES directly reflect mean differences in milliseconds.

This table corresponds to Table S2 in the supplemental material.

# Run model
LMM_pes <- lmer(pes ~ stimulation * group +
  (1 + stimulation | participant_id),
data = single_trial_data_clean[!is.na(single_trial_data_clean$pes), ],
REML = TRUE,
control = lmerControl(optimizer = "bobyqa")
)


# Check model output and PCA of random-effects variance-covariance estimates
# summary(LMM_pes)        # Model does converge
# isSingular(LMM_pes)     # No singular fit
# summary(rePCA(LMM_pes)) # All terms explain variance


# Display results (fixed effects)
tab_model(LMM_pes,
  dv.labels = "PES", pred.labels = labels, show.stat = TRUE, show.icc = FALSE, show.r2 = FALSE,
  show.re.var = FALSE, show.ngroups = FALSE, minus.sign = "−", string.pred = "Fixed effects", 
  string.est = "b", string.stat = "t", string.ci = "95 % CI", p.val = "satterthwaite", 
  title = "Results of the LMM Predicting PES as a Function of Stimulation Condition 
  (Cathodal - Sham) and Group (OCD - Healthy Controls)"
)


# Display random effects
print("Random effects:")
print(VarCorr(LMM_pes), digits = 4, comp = "Std.Dev.")

Results of the LMM Predicting PES as a Function of Stimulation Condition (Cathodal - Sham) and Group (OCD - Healthy Controls)

	PES
Fixed effects	b	95 % CI	t	p
Intercept	35.71	29.40 – 42.02	11.09	<0.001
Stimulation	1.08	-8.73 – 10.89	0.22	0.830
Group	9.59	-3.03 – 22.21	1.49	0.142
Stimulation x Group	-3.08	-22.70 – 16.54	-0.31	0.760
Observations	2622

[1] "Random effects:"
 Groups         Name           Std.Dev. Corr 
 participant_id (Intercept)    19.17         
                stimulation2-1 23.80    -0.02
 Residual                      95.82         

The LMM on PES yielded a significant intercept, reflecting that participants slowed down after error commission. The groups did not differ significantly in PES. There was no significant main effect of tDCS condition and no interaction between group and tDCS condition.

Response Accuracy

Response accuracy was analyzed using a binomial generalized linear mixed model (GLMM). For the GLMM on accuracy, estimates reflect odds ratios for a correct response.

This table corresponds to Table S2 in the supplemental material.

# Run model
GLMM_accuracy <- glmer(accuracy_numeric ~ stimulation * group +
  (1 + stimulation | participant_id),
data = single_trial_data_clean,
family = binomial,
control = glmerControl(optimizer = "bobyqa")
)


# Check model output and PCA of random-effects variance-covariance estimates
# summary(GLMM_accuracy)        # Model does converge
# isSingular(GLMM_accuracy)     # No singular fit
# summary(rePCA(GLMM_accuracy)) # All terms explain variance


# Display results (fixed effects)
tab_model(GLMM_accuracy,
  dv.labels = "Accuracy", pred.labels = labels, show.stat = TRUE, show.icc = FALSE, show.r2 = FALSE,
  show.re.var = FALSE, show.ngroups = FALSE, minus.sign = "−", string.pred = "Fixed effects", 
  string.est = "Odds ratio", string.stat = "z", string.ci = "95 % CI", title = "Results of the Binomial GLMM 
  Predicting Response Accuracy as a Function of Stimulation Condition (Cathodal - Sham) and 
  Group (OCD - Healthy Controls)"
)


# Display random effects
print("Random effects:")
print(VarCorr(GLMM_accuracy), digits = 2, comp = "Std.Dev.")

Results of the Binomial GLMM Predicting Response Accuracy as a Function of Stimulation Condition (Cathodal - Sham) and Group (OCD - Healthy Controls)

	Accuracy
Fixed effects	Odds ratio	95 % CI	z	p
Intercept	16.45	14.84 – 18.24	53.10	<0.001
Stimulation	0.97	0.88 – 1.07	-0.56	0.573
Group	1.18	0.96 – 1.44	1.54	0.123
Stimulation x Group	0.84	0.70 – 1.02	-1.71	0.087
Observations	53093

[1] "Random effects:"
 Groups         Name           Std.Dev. Corr
 participant_id (Intercept)    0.37         
                stimulation2-1 0.24     0.01

Results of the GLMM indicated that there was no significant difference in response accuracy between groups or tDCS conditions and no significant interaction between group and tDCS condition.

References

---

Bates, D., Kliegl, R., Vasishth, S., & Baayen, H. (2015). Parsimonious mixed models. arXiv. https://arxiv.org/abs/1506.04967v2

Box, G. E., & Cox, D. R. (1964). An analysis of transformations. Journal of the Royal Statistical Society: Series B (Methodological), 26(2), 211-243. https://doi.org/10.1111/j.2517-6161.1964.tb00553.x

Dutilh, G., van Ravenzwaaij, D., Nieuwenhuis, S., van der Maas, H. L. J., Forstmann, B. U., & Wagenmakers, E.-J. (2012). How to measure post-error slowing: A confound and a simple solution. Journal of Mathematical Psychology, 56(3), 208-216. https://doi.org/10.1016/j.jmp.2012.04.001

Morey, R. (2008). Confidence intervals from normalized data: A correction to Cousineau (2005). Tutorials in Quantitative Methods for Psychology, 4(2), 61-64. https://doi.org/10.20982/tqmp.04.2.p061

Session Info

---

sessionInfo()

R version 3.6.1 (2019-07-05)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19044)

Matrix products: default

locale:
[1] LC_COLLATE=German_Germany.1252  LC_CTYPE=German_Germany.1252   
[3] LC_MONETARY=German_Germany.1252 LC_NUMERIC=C                   
[5] LC_TIME=German_Germany.1252    

attached base packages:
[1] stats     graphics  grDevices datasets  utils     methods   base     

other attached packages:
[1] plyr_1.8.6       sjPlot_2.8.6     lmerTest_3.1-3   lme4_1.1-25     
[5] Matrix_1.2-17    MASS_7.3-53      kableExtra_1.3.1 knitr_1.30      
[9] dplyr_1.0.6     

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.5          mvtnorm_1.1-1       lattice_0.20-38    
 [4] tidyr_1.1.2         digest_0.6.27       utf8_1.1.4         
 [7] R6_2.5.0            backports_1.2.0     evaluate_0.14      
[10] highr_0.8           httr_1.4.2          ggplot2_3.3.2      
[13] pillar_1.6.0        rlang_0.4.11        rstudioapi_0.11    
[16] minqa_1.2.4         performance_0.7.1.1 nloptr_1.2.2.2     
[19] effectsize_0.4.0    rmarkdown_2.5       ggeffects_0.16.0   
[22] splines_3.6.1       webshot_0.5.2       statmod_1.4.35     
[25] stringr_1.4.0       munsell_0.5.0       broom_0.7.2        
[28] modelr_0.1.8        compiler_3.6.1      numDeriv_2016.8-1.1
[31] xfun_0.19           pkgconfig_2.0.3     parameters_0.13.0  
[34] htmltools_0.5.0     insight_0.14.0      tidyselect_1.1.0   
[37] tibble_3.0.4        codetools_0.2-16    fansi_0.4.1        
[40] viridisLite_0.3.0   crayon_1.3.4        sjmisc_2.8.5       
[43] grid_3.6.1          nlme_3.1-140        xtable_1.8-4       
[46] gtable_0.3.0        lifecycle_1.0.0     magrittr_2.0.1     
[49] bayestestR_0.9.0    scales_1.1.1        estimability_1.3   
[52] stringi_1.5.3       renv_0.12.0         xml2_1.3.2         
[55] ellipsis_0.3.2      generics_0.1.0      vctrs_0.3.8        
[58] boot_1.3-22         sjlabelled_1.1.7    tools_3.6.1        
[61] glue_1.4.2          purrr_0.3.4         sjstats_0.18.0     
[64] emmeans_1.5.2-1     yaml_2.2.1          colorspace_1.4-1   
[67] rvest_1.0.2