minor edits and rearranging

code/highspeed-analysis-behavior.Rmd

@@ -1,50 +1,11 @@
----
-title: "Faster than thought: Detecting sub-second activation sequences with sequential fMRI pattern analysis"
-subtitle: "Behavioral results"
-author: "Lennart Wittkuhn & Nicolas W. Schuck"
-date: "Last update: `r format(Sys.time(), '%d %B, %Y')`"
-site: bookdown::bookdown_site
-output:
-  bookdown::gitbook:
-    config:
-      toc:
-        collapse: subsection
-        scroll_highlight: yes
-        before: null
-        after: null
-      toolbar:
-        position: fixed
-      edit : null
-      download: null
-      search: yes
-      fontsettings:
-        theme: white
-        family: sans
-        size: 2
-      sharing:
-        facebook: yes
-        github: yes
-        twitter: yes
-        linkedin: no
-        weibo: no
-        instapaper: no
-        vk: no
-        all: ['facebook', 'twitter', 'linkedin']
-      info: yes
-      code_folding: "hide"
-fig.align: "center"
-header-includes:
-  - \usepackage{fontspec}
-  - \setmainfont{AgfaRotisSansSerif}
-email: wittkuhn@mpib-berlin.mpg.de
-documentclass: book
-link-citations: yes
----
-
 # Behavior
-## Initialization {.tabset .tabset-fade .tabset-pills}
+
+## Initialization
 
 ### Load data and files
+
+We set the paths and source the basic setup script:
+
 ```{r, warning=FALSE, message=FALSE}
 # find the path to the root of this project:
 if (!requireNamespace("here")) install.packages("here")
@@ -57,7 +18,10 @@ if ( basename(here::here()) == "highspeed" ) {
 source(file.path(path_root, "code", "highspeed-analysis-setup.R"))
 ```
 
-### Assign signal-detection labels
+### Signal-detection labeling
+
+We assign labels from signal detection theory that will be used in one of the analyses below:
+
 ```{r}
 # denotes misses (key was not pressed and stimulus was upside-down):
 dt_events$sdt_type[
@@ -74,7 +38,9 @@ dt_events$sdt_type[
 ```
 
 ## Stimulus timings
+
 We calculate the differences between consecutive stimulus onsets:
+
 ```{r}
 dt_events %>%
   # get duration of stimuli by calculating differences between consecutive onsets:
@@ -159,9 +125,12 @@ dt_odd_iti_mean = dt_events %>%
 rmarkdown::paged_table(dt_odd_iti_mean)
 ```
 
-## Behavioral performance
+## Overview: Behavioral performance
+
+### Mean accuracy
+
+We calculate the mean behavioral accuracy across all trials of all three task conditions (slow, sequence, and repetition trials):
 
-### Overview: Mean accuracy
 ```{r, echo = TRUE}
 chance_level = 50
 dt_acc = dt_events %>%
@@ -193,32 +162,42 @@ dt_acc = dt_events %>%
 rmarkdown::paged_table(dt_acc)
 ```
 
-```{r, echo = TRUE}
+We create a list of participants that will be excluded because their performance is below the 50% chance level in either or both sequence and repetition trials:
+
+```{r, echo=TRUE}
 # create a list with all excluded subject ids and print the list:
 subjects_excluded = unique(dt_acc$subject[dt_acc$exclude == "yes"])
 print(subjects_excluded)
 ```
 
-```{r, echo = TRUE, results = "hold"}
+We calculate the mean behavioral accuracy across all three task conditions (slow, sequence, and repetition trials), *excluding* partipants that performed below chance on either or both sequence and repetition trials:
+
+```{r, echo=TRUE, results="hold"}
 dt_acc_mean = dt_acc %>%
   # filter out all data of excluded participants:
   filter(!(subject %in% unique(subject[exclude == "yes"]))) %>%
   # check if the number of participants matches expectations:
-  verify(length(unique(subject)) == 36) %>% setDT(.) %>%
+  verify(length(unique(subject)) == 36) %>%
+  setDT(.) %>%
   # calculate mean behavioral accuracy across participants for each condition:
   .[, by = .(condition), {
-    ttest_results = t.test(mean_accuracy, mu = chance_level, alternative = "greater")
+    ttest_results = t.test(
+      mean_accuracy, mu = chance_level, alternative = "greater")
     list(
-      pvalue = round_pvalues(ttest_results$p.value),
+      pvalue = ttest_results$p.value,
+      pvalue_rounded = round_pvalues(ttest_results$p.value),
       tvalue = round(ttest_results$statistic, digits = 2),
+      conf_lb = round(ttest_results$conf.int[1], digits = 2),
+      conf_ub = round(ttest_results$conf.int[2], digits = 2),
       df = ttest_results$parameter,
       num_subs = .N,
-      mean_accuracy = mean(mean_accuracy),
-      SD = sd(mean_accuracy),
+      mean_accuracy = round(mean(mean_accuracy), digits = 2),
+      SD = round(sd(mean_accuracy), digits = 2),
       cohens_d = round((mean(mean_accuracy) - chance_level) / sd(mean_accuracy), 2),
       sem_upper = mean(mean_accuracy) + (sd(mean_accuracy)/sqrt(.N)),
       sem_lower = mean(mean_accuracy) - (sd(mean_accuracy)/sqrt(.N))
-    )}] %>% verify(num_subs == 36) %>%
+    )}] %>%
+  verify(num_subs == 36) %>%
   # create a short name for the conditions:
   mutate(condition_short = substr(condition, start = 1, stop = 3)) %>%
   # reorder the condition factor in descending order of accuracy:
@@ -228,7 +207,9 @@ rmarkdown::paged_table(dt_acc_mean)
 ```
 
 ### Above-chance performance
+
 We plot only data of above-chance performers:
+
 ```{r, echo = FALSE}
 fig_behav_all = ggplot(data = subset(dt_acc, exclude == "no"), aes(
   x = as.factor(condition_short), y = as.numeric(mean_accuracy),
@@ -249,7 +230,9 @@ fig_behav_all
 ```
 
 ### Below chance performance
+
 We plot data of all participants with below chance performers highlighted in red.
+
 ```{r, echo = FALSE}
 fig_behav_all_outlier = ggplot(data = dt_acc_mean,
   mapping = aes(x = as.factor(condition_short), y = as.numeric(mean_accuracy),
@@ -279,29 +262,35 @@ fig_behav_all_outlier = ggplot(data = dt_acc_mean,
 fig_behav_all_outlier
 ```
 
-## Oddball task
-### Mean accuracy
-Accuracy on oddball trials across all trials (in final sample):
-```{r}
+## Slow trials
+
+### Mean accuracy (all trials)
+
+We calculate the mean accuracy on slow trials (oddball task condition) across all trials in the final sample (only participants who performed above chance):
+
+```{r, echo=TRUE}
+# we use the dataframe containing the accuracy data
 dt_acc_odd = dt_acc %>%
   # filter for oddball / slow trials only:
   filter(condition == "oddball") %>%
   # exclude participants with below chance performance::
   filter(!(subject %in% subjects_excluded)) %>%
-  # verify that the number of participants is correct:
+  # verify that the number of participants (final sample) is correct:
   verify(all(.N == 36))
 ```
 
-```{r, echo = FALSE, fig.width=3}
+We plot the mean behavioral accuracy on slow trials (oddball task condition) in the final sample:
+
+```{r, echo=TRUE, fig.width=3}
 fig_behav_odd = ggplot(data = dt_acc_odd, aes(
   x = "mean_acc", y = as.numeric(mean_accuracy))) +
   geom_bar(stat = "summary", fun = "mean", fill = "lightgray") +
   #geom_dotplot(binaxis = "y", stackdir = "center", stackratio = 0.5,
   #             color = "black", fill = "lightgray", alpha = 0.5,
   #             inherit.aes = TRUE, binwidth = 0.5) +
-  #geom_point(position = position_jitter(width = 0.2, height = 0, seed = 2),
-  #           alpha = 0.5, inherit.aes = TRUE, pch = 21,
-  #           color = "black", fill = "lightgray") +
+  geom_point(position = position_jitter(width = 0.2, height = 0, seed = 2),
+             alpha = 0.5, inherit.aes = TRUE, pch = 21,
+             color = "black", fill = "lightgray") +
   geom_errorbar(stat = "summary", fun.data = "mean_se", width = 0.0, color = "black") +
   ylab("Accuracy (%)") + xlab("Condition") +
   scale_color_manual(values = c("darkgray", "red"), name = "Outlier") +
@@ -315,7 +304,10 @@ fig_behav_odd = ggplot(data = dt_acc_odd, aes(
 fig_behav_odd
 ```
 
-### Accuracy across runs
+### Mean accuracy (per run)
+
+We calculate the mean behavioral accuracy on slow trials (oddball task condition) for each of the eight task runs *for each* participant:
+
 ```{r}
 # calculate the mean accuracy per session and run for every participant:
 dt_odd_behav_run_sub = dt_events %>%
@@ -331,11 +323,13 @@ dt_odd_behav_run_sub = dt_events %>%
     mean_accuracy = mean(accuracy))] %>%
   # express accuracy in percent by multiplying with 100:
   transform(mean_accuracy = mean_accuracy * 100) %>%
-  # z-score the accuracy values:
-  #mutate(mean_accuracy_z = scale(mean_accuracy, scale = TRUE, center = TRUE)) %>%
   # check whether the mean accuracy is within the expected range of 0 to 100:
   assert(within_bounds(lower.bound = 0, upper.bound = 100), mean_accuracy)
+```
+
+We calculate the mean behavioral accuracy on slow trials (oddball task condition) for each of the eight task runs *across* participants:
 
+```{r}
 # calculate mean accuracy per session and run across participants:
 dt_odd_behav_run_mean = dt_odd_behav_run_sub %>%
   setDT(.) %>%
@@ -345,14 +339,15 @@ dt_odd_behav_run_mean = dt_odd_behav_run_sub %>%
     num_subs = .N,
     sem_upper = mean(mean_accuracy) + (sd(mean_accuracy)/sqrt(.N)),
     sem_lower = mean(mean_accuracy) - (sd(mean_accuracy)/sqrt(.N))
-  )] %>% verify(num_subs == 36) %>%
+  )] %>%
+  verify(num_subs == 36) %>%
   # z-score the accuracy values:
   mutate(mean_accuracy_z = scale(mean_accuracy, scale = TRUE, center = TRUE))
 ```
 
+We run a LME model to test the linear effect of task run on behavioral accuracy:
 
-We run a LME model to test the linear effect of experiment run on accuracy:
-```{r, results = "hold"}
+```{r, results="hold", echo=TRUE}
 lme_odd_behav_run = lmerTest::lmer(
   mean_accuracy ~ run_study + (1 + run_study | subject),
   data = dt_odd_behav_run_sub, na.action = na.omit, control = lcctrl)
@@ -361,10 +356,14 @@ anova(lme_odd_behav_run)
 ```
 
 We run a second model to test run- and session-specific effects:
-```{r, results = "hold"}
+
+```{r, results="hold"}
 dt <- dt_odd_behav_run_sub %>%
   transform(run_session = as.factor(paste0("run-0", run_session)),
             session = as.factor(paste0("ses-0", session)))
+```
+
+```{r, results="hold"}
 lme_odd_behav_run = lmerTest::lmer(
   mean_accuracy ~ session + run_session + (1 + session + run_session | subject),
   data = dt, na.action = na.omit, control = lcctrl)
@@ -374,7 +373,9 @@ anova(lme_odd_behav_run)
 rm(dt)
 ```
 
-```{r, echo=FALSE}
+We plot the behavioral accuracy on slow trials (oddball task condition) across task runs (x-axis) for each study session (panels):
+
+```{r, echo=TRUE}
 # change labels of the facet:
 facet_labels_new = unique(paste0("Session ", dt_events$session))
 facet_labels_old = as.character(unique(dt_events$session))
@@ -393,9 +394,11 @@ plot_odd_run = ggplot(data = dt_odd_behav_run_mean, mapping = aes(
 plot_odd_run
 ```
 
-### Performance for misses and false alarms
+### Misses vs. false alarms
 
-```{r}
+We calculate the mean frequency of misses (missed response to upside-down images) and false alarms (incorrect response to upright images):
+
+```{r, echo=TRUE}
 dt_odd_behav_sdt_sub = dt_events %>%
   # exclude participants performing below chance:
   filter(!(subject %in% subjects_excluded)) %>%
@@ -419,6 +422,8 @@ dt_odd_behav_sdt_sub = dt_events %>%
   mutate(sdt_type_numeric = ifelse(sdt_type == "false alarm", 1, -1))
 ```
 
+We run a LME model to test the effect of signal detection type (miss vs. false alarm), task run and session on the frequency of those events:
+
 ```{r, results = "hold"}
 lme_odd_behav_sdt = lmer(
   freq ~ sdt_type + run_session * session + (1 + run_session + session | subject),
@@ -430,7 +435,9 @@ emmeans_pvalues = round_pvalues(summary(emmeans_results[[2]])$p.value)
 emmeans_results
 ```
 
-```{r, echo = FALSE}
+We plot the frequency of misses and false alarms as a function of task run and study session:
+
+```{r, echo=TRUE}
 plot_odd_sdt = ggplot(data = dt_odd_behav_sdt_sub, mapping = aes(
   y = as.numeric(freq), x = as.numeric(run_session),
   fill = as.factor(sdt_type))) +
@@ -438,7 +445,7 @@ plot_odd_sdt = ggplot(data = dt_odd_behav_sdt_sub, mapping = aes(
   stat_summary(geom = "errorbar", fun.data = mean_se, position = position_dodge(0.9), width = 0) +
   facet_wrap(~ as.factor(session), labeller = as_labeller(facet_labels_new)) +
   ylab("Frequency (%)") + xlab("Run") +
-  ylim(c(0, 10)) +
+  #ylim(c(0, 10)) +
   #coord_capped_cart(left = "both", bottom = "both", expand = TRUE, ylim = c(0,10)) +
   scale_fill_viridis(name = "Error", discrete = TRUE) +
   theme(axis.ticks.x = element_blank(), axis.line.x = element_blank()) +
@@ -451,9 +458,11 @@ plot_odd_sdt
 
 ## Sequence trials
 
-### Mean accuracy
+### Effect of sequence speed
 
-```{r}
+We calculate the mean behavioral accuracy on sequence trials for each of the five sequence speeds (inter-stimulus intervals):
+
+```{r, echo=TRUE}
 dt_seq_behav = dt_events %>%
   # filter behavioral events data for sequence trials only:
   filter(condition == "sequence") %>%
@@ -476,8 +485,6 @@ dt_seq_behav = dt_events %>%
   setDT(.)
 ```
 
-### Effect of sequence speed
-
 ```{r}
 dt_seq_behav_speed = dt_seq_behav %>%
   # filter out excluded subjects:
@@ -499,7 +506,9 @@ dt_seq_behav_speed = dt_seq_behav %>%
   setDT(.)
 ```
 
-```{r, results="hold"}
+We run a LME model to test the effect of sequence speed (inter-stimulus interval) on mean behavioral accuracy on sequence trials:
+
+```{r, results="hold", echo=TRUE}
 lme_seq_behav = lmer(
   mean_accuracy ~ trial_speed + (1 + trial_speed | subject),
   data = dt_seq_behav_speed, na.action = na.omit, control = lcctrl)
@@ -511,7 +520,7 @@ emmeans_results
 emmeans_pvalues
 ```
 
-### Difference from chance
+We compare mean behavioral accuracy at each sequence speed level to the chancel level of 50%:
 
 ```{r}
 chance_level = 50
@@ -523,14 +532,17 @@ dt_seq_behav_mean = dt_seq_behav_speed %>%
   list(
     mean_accuracy = round(mean(mean_accuracy), digits = 2),
     sd_accuracy = round(sd(mean_accuracy), digits = 2),
-    tvalue = round(ttest_results$estimate, digits = 2),
+    tvalue = round(ttest_results$statistic, digits = 2),
+    conf_lb = round(ttest_results$conf.int[1], digits = 2),
+    conf_ub = round(ttest_results$conf.int[2], digits = 2),
     pvalue = ttest_results$p.value,
     cohens_d = round((mean(mean_accuracy) - chance_level)/sd(mean_accuracy), 2),
     df = ttest_results$parameter,
     num_subs = .N,
     sem_upper = mean(mean_accuracy) + (sd(mean_accuracy)/sqrt(.N)),
     sem_lower = mean(mean_accuracy) - (sd(mean_accuracy)/sqrt(.N))
-  )}] %>% verify(num_subs == 36) %>%
+  )}] %>%
+  verify(num_subs == 36) %>%
   mutate(sem_range = sem_upper - sem_lower) %>%
   mutate(pvalue_adjust = p.adjust(pvalue, method = "fdr")) %>%
   mutate(pvalue_adjust_round = round_pvalues(pvalue_adjust))
@@ -538,7 +550,7 @@ dt_seq_behav_mean = dt_seq_behav_speed %>%
 rmarkdown::paged_table(dt_seq_behav_mean)
 ```
 
-### Reduction in accuracy
+We calculate the reduction in mean behavioral accuracy comparing the fastest and slowest speed condition:
 
 ```{r}
 a = dt_seq_behav_mean$mean_accuracy[dt_seq_behav_mean$trial_speed == 2.048]
@@ -552,9 +564,9 @@ fig_seq_speed = ggplot(data = dt_seq_behav_speed, mapping = aes(
   y = as.numeric(mean_accuracy), x = as.factor(as.numeric(trial_speed)*1000),
   fill = as.factor(trial_speed), color = as.factor(trial_speed))) +
   geom_bar(stat = "summary", fun = "mean") +
-  #geom_dotplot(binaxis = "y", stackdir = "center", stackratio = 0.5,
-  #             color = "black", alpha = 0.5,
-  #             inherit.aes = TRUE, binwidth = 2) +
+  geom_dotplot(binaxis = "y", stackdir = "center", stackratio = 0.5,
+               color = "black", alpha = 0.5,
+               inherit.aes = TRUE, binwidth = 2) +
   #geom_point(position = position_jitter(width = 0.2, height = 0, seed = 3),
   #           alpha = 0.5, pch = 21, color = "black") +
   geom_errorbar(stat = "summary", fun.data = "mean_se", width = 0.0, color = "black") +
@@ -601,9 +613,9 @@ fig_seq_position = ggplot(data = dt_seq_behav_position, mapping = aes(
   y = as.numeric(mean_accuracy), x = as.factor(trial_target_position),
   fill = as.factor(trial_target_position), color = as.factor(trial_target_position))) +
   geom_bar(stat = "summary", fun = "mean") +
-  #geom_dotplot(binaxis = "y", stackdir = "center", stackratio = 0.5,
-  #             color = "black", alpha = 0.5,
-  #             inherit.aes = TRUE, binwidth = 2) +
+  geom_dotplot(binaxis = "y", stackdir = "center", stackratio = 0.5,
+               color = "black", alpha = 0.5,
+               inherit.aes = TRUE, binwidth = 2) +
   #geom_point(pch = 21, alpha = 0.5, color = "black",
   #           position = position_jitter(height = 0, seed = 3)) +
   geom_errorbar(stat = "summary", fun.data = "mean_se", width = 0.0, color = "black") +
@@ -620,6 +632,10 @@ fig_seq_position
 
 ## Repetition trials
 
+### Mean accuracy
+
+We calculate mean behavioral accuracy in repetition trials for each participant:
+
 ```{r}
 dt_rep_behav = dt_events %>%
   # filter for repetition trials only:
@@ -642,10 +658,14 @@ dt_rep_behav = dt_events %>%
   # transform mean accuracy into percent (%)
   transform(mean_accuracy = mean_accuracy * 100) %>%
   # check if accuracy values range between 0 and 100
-  verify(between(x = mean_accuracy, lower = 0, upper = 100))
+  verify(between(x = mean_accuracy, lower = 0, upper = 100)) %>%
+  mutate(interference = ifelse(
+    trial_target_position == 2, "fwd", trial_target_position)) %>%
+  transform(interference = ifelse(
+    trial_target_position == 9, "bwd", interference))
 ```
 
-### Difference from chance
+We run separate one-sided one-sample t-tests to access if mean behavioral performance for every repetition condition differs from a 50% chance level:
 
 ```{r}
 chance_level = 50
@@ -660,8 +680,10 @@ dt_rep_behav_chance = dt_rep_behav %>%
   list(
     mean_accuracy = round(mean(mean_accuracy), digits = 2),
     sd_accuracy = round(sd(mean_accuracy), digits = 2),
-    tvalue = round(ttest_results$estimate, digits = 2),
+    tvalue = round(ttest_results$statistic, digits = 2),
     pvalue = ttest_results$p.value,
+    conf_lb = round(ttest_results$conf.int[1], digits = 2),
+    conf_ub = round(ttest_results$conf.int[2], digits = 2),
     cohens_d = round((mean(mean_accuracy) - chance_level)/sd(mean_accuracy), 2),
     df = ttest_results$parameter,
     num_subs = .N,
@@ -695,11 +717,14 @@ fig_behav_rep = ggplot(data = plot_data, mapping = aes(
   y = as.numeric(mean_accuracy), x = fct_rev(as.factor(interference)),
   fill = as.factor(interference))) +
   geom_bar(stat = "summary", fun = "mean") +
-  #geom_dotplot(data = dt_rep_behav %>%
-  #               filter(!(subject %in% subjects_excluded)) %>%
-  #               filter(trial_target_position %in% c(2,9)),
-  #             binaxis = "y", stackdir = "center", stackratio = 0.5,
-  #             color = "black", alpha = 0.5, inherit.aes = TRUE, binwidth = 2) +
+  geom_dotplot(data = dt_rep_behav %>%
+                 filter(!(subject %in% subjects_excluded)) %>%
+                 filter(trial_target_position %in% c(2,9)) %>% setDT(.) %>%
+                 .[, by = .(subject, interference), .(
+                   mean_accuracy = mean(mean_accuracy)
+                 )],
+               binaxis = "y", stackdir = "center", stackratio = 0.5,
+               color = "black", alpha = 0.5, inherit.aes = TRUE, binwidth = 2) +
   geom_errorbar(aes(ymin = sem_lower, ymax = sem_upper), width = 0.0, color = "black") +
   ylab("Accuracy (%)") + xlab("Interfererence") +
   ggtitle("Repetition") +
@@ -720,6 +745,17 @@ plot_behav_rep_all = ggplot(data = plot_data, mapping = aes(
   y = as.numeric(mean_accuracy), x = as.numeric(trial_target_position),
   fill = as.numeric(trial_target_position))) +
   geom_bar(stat = "identity") +
+  geom_dotplot(data = dt_rep_behav %>%
+                 filter(!(subject %in% subjects_excluded)) %>%
+                 filter(trial_target_position %in% seq(2,9)) %>% setDT(.) %>%
+                 .[, by = .(subject, trial_target_position), .(
+                   mean_accuracy = mean(mean_accuracy)
+                 )],
+               aes(x = as.numeric(trial_target_position),
+                   fill = as.numeric(trial_target_position),
+                   group = as.factor(trial_target_position)),
+               binaxis = "y", stackdir = "center", stackratio = 0.5,
+               inherit.aes = TRUE, binwidth = 2, color = "black", alpha = 0.5) +
   geom_errorbar(aes(ymin = sem_lower, ymax = sem_upper), width = 0.0, color = "black") +
   geom_text(aes(y = as.numeric(mean_accuracy) + 10, label = cohens_d), size = 2.5) +
   ylab("Accuracy (%)") + xlab("First / second item repetitions") +