Preliminaries.

rm(list=ls())
suppressPackageStartupMessages(library(dplyr))
suppressPackageStartupMessages(library(langcog))
suppressPackageStartupMessages(library(lme4))
suppressPackageStartupMessages(library(lmerTest))
library(ggplot2)
library(magrittr)
library(knitr)
library(tidyr)
## 
## Attaching package: 'tidyr'
## The following object is masked from 'package:magrittr':
## 
##     extract
## The following object is masked from 'package:Matrix':
## 
##     expand
opts_chunk$set(fig.width=8, fig.height=5, 
                      echo=TRUE, warning=FALSE, message=FALSE, cache=TRUE)
theme_set(theme_bw())

1 Experiment 1: Experts with dual task

1.1 Data prep

Read data. Coding info from Katie.

  • Heaven/Earth 0=Heaven; 1=Earth
  • Condition: 1= In Play; 2=Out of play; 3=leading; 4=trailing
  • Number Requested,
  • X: 1=Left most column, counting to the right to max of 2 or 3“,
  • Y: 1 equals top row (i.e., out of play heavenly bead), counting downwards to 7.
  • Correct on Search task? Trials weith incorrect search responses were excluded from analysis
  • RT
  • Abacus entered: Only relevent to dual task,
  • Correct on Abacus: only relevent to dual task
  • “Number of columns: occasionally you will see 0s in this column, this means the participant was in a pilot version with 4 columns and should be excluded.”
  • Subject Number: Assigned when extracted from .mat instead of alpha numeric number to make some things in life easier,
  • <3Std Dev: removing outliers (2s). Calculated in linear space. Needs to be redone in log space,
  • Inout or Leading trailing trial,
  • Expertise level: ranges from 0 (none) to 2(has used abacus)
dual_experts_raw <- read.csv("data/Upright dual data experts w. subject IDs and order.csv")

names(dual_experts_raw) <- c("bead_type", "condition", "number_requested", 
                             "X_pos","Y_pos","search_correct","RT",
                             "abacus_val","abacus_correct","n_col",
                             "subnum","outlier","trial_type",
                             "subid","n_conditions","order","subid_xexpts")
dual_experts_raw %<>% 
  mutate(bead_type = factor(bead_type, 
                            levels = c(0,1), 
                            labels = c("heaven","earth")),
         condition = factor(condition, 
                            levels = c(1,2,3,4), 
                            labels = c("in play", "out of play", 
                                       "leading","trailing")))

Exclusions. Filter pilot participants.

pilot_subs <- dual_experts_raw %>%
  group_by(subnum) %>%
  summarise(pilot = any(n_col == 0)) %>%
  filter(pilot) 
  
dual_experts <- filter(dual_experts_raw, 
                       !subnum %in% pilot_subs$subnum)
dual_experts %<>% 
  group_by(subnum) %>%
  mutate(trial_num = 1:n())

Check to make sure we have a consistent number of trials, no training trials.

qplot(subnum, trial_num, data = dual_experts)

In this dataset we appear to be missing the end of a few participants.

Next, RT exclusions. Note that there are a few 0 RTs. What’s the deal with these?

sum(dual_experts$RT == 0)
## [1] 37
dual_experts %<>% filter(RT > 0, 
                         !is.na(RT))

Linear space.

qplot(RT, data = dual_experts, 
      fill = RT > mean(RT) + 3*sd(RT))

mean(dual_experts$RT)
## [1] 3.860437
median(dual_experts$RT)
## [1] 3.31

Log space looks better.

qplot(log(RT), data = dual_experts, 
      fill = log(RT) > mean(log(RT)) + 3*sd(log(RT)) |
        log(RT) < mean(log(RT)) - 3*sd(log(RT)))

Clip these.

lmean <- mean(log(dual_experts$RT))
lsd <- sd(log(dual_experts$RT))
dual_experts$RT[log(dual_experts$RT) > lmean + 3*lsd |
                  log(dual_experts$RT) < lmean - 3*lsd] <- NA

Replot in linear space just to check.

qplot(RT, data = dual_experts)

Looks good.

1.2 Demographics and descriptives

kable(dual_experts %>% 
        group_by(subnum) %>% 
        summarise(n = n()) %>% 
        summarise(n = n()))

1.3 n

63

kable(dual_experts %>% 
        group_by(subnum) %>% 
        summarise(RT = mean(RT, na.rm=TRUE)) %>% 
        mutate(condition = "all") %>% group_by(condition) %>%
        multi_boot_standard(col = "RT"))
condition mean ci_lower ci_upper
all 3.756448 3.548135 3.984962

1.4 RT and accuracy analyses

Basic analyses.

ms <- dual_experts %>%
  filter(abacus_correct == 1, 
         search_correct == 1) %>%
  group_by(subnum, condition) %>%
  summarise(RT = mean(RT, na.rm=TRUE)) %>%
  group_by(condition) %>%
  multi_boot_standard(col = "RT")

kable(ms, digits = 2)
condition mean ci_lower ci_upper
in play 3.63 3.40 3.86
out of play 3.88 3.66 4.10
leading 3.86 3.63 4.09
trailing 3.66 3.43 3.89 
ggplot(ms,aes(x = condition, y = mean, fill = condition)) + 
  geom_bar(stat = "identity") + 
  geom_linerange(aes(ymin = ci_lower, ymax = ci_upper))

Add two other variables: bead type and number of columns.

ms <- dual_experts %>%
  filter(abacus_correct == 1, 
         search_correct == 1) %>%
  group_by(subnum, condition, bead_type, n_col) %>%
  summarise(RT = mean(RT, na.rm=TRUE)) %>%
  group_by(condition, bead_type, n_col) %>%
  multi_boot_standard(col = "RT")

kable(ms, digits = 2)
condition bead_type n_col mean ci_lower ci_upper
in play heaven 2 3.39 3.13 3.69
in play heaven 3 4.40 4.00 4.90
in play earth 2 3.07 2.87 3.31
in play earth 3 4.03 3.75 4.32
out of play heaven 2 3.74 3.45 4.06
out of play heaven 3 4.46 4.15 4.80
out of play earth 2 3.44 3.24 3.65
out of play earth 3 4.18 3.91 4.44
leading heaven 2 3.97 3.66 4.29
leading heaven 3 4.41 4.12 4.71
leading earth 2 3.51 3.29 3.75
leading earth 3 4.00 3.77 4.27
trailing heaven 2 3.71 3.43 4.04
trailing heaven 3 4.36 4.04 4.69
trailing earth 2 3.16 2.98 3.33
trailing earth 3 3.91 3.64 4.21 
ggplot(ms,aes(x = condition, y = mean, fill = condition)) + 
  geom_bar(stat = "identity") + 
  geom_linerange(aes(ymin = ci_lower, ymax = ci_upper)) + 
  facet_grid(bead_type ~ n_col)

It’s clear that the effects are being driven by the 2-column displays, and especiallywith the earthly beads. (Though there are probably fewer heavenly bead trials, no?).

Order effects.

ms <- dual_experts %>%
  filter(abacus_correct == 1, 
         search_correct == 1) %>%
  mutate(second_half = trial_num > 64) %>%
  group_by(subnum, condition, second_half) %>%
  summarise(RT = mean(RT, na.rm=TRUE)) %>%
  group_by(condition, second_half) %>%
  multi_boot_standard(col = "RT")

kable(ms, digits = 2)
condition second_half mean ci_lower ci_upper
in play FALSE 3.88 3.60 4.14
in play TRUE 3.37 3.12 3.61
out of play FALSE 4.12 3.88 4.38
out of play TRUE 3.69 3.46 3.95
leading FALSE 4.14 3.92 4.39
leading TRUE 3.59 3.36 3.83
trailing FALSE 3.97 3.72 4.24
trailing TRUE 3.35 3.15 3.57 
ggplot(ms,aes(x = condition, y = mean, fill = condition)) + 
  geom_bar(stat = "identity") + 
  facet_wrap(~second_half) + 
  geom_linerange(aes(ymin = ci_lower, ymax = ci_upper))

1.5 Stats

Basic LMER confirms highly significant effects for in/out of play, and leading/trailing. Could do better parameterization, but this is still very clear. More random effects don’t converge.

kable(summary(lmer(log(RT) ~ trial_num + condition + 
                     (condition | subnum), 
                   data = filter(dual_experts, 
                                 search_correct == 1, 
                                 abacus_correct == 1)))$coefficients, digits = 3)
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 1.330 0.035 68.048 38.550 0.000
trial_num -0.003 0.000 7345.036 -22.050 0.000
conditionout of play 0.092 0.015 64.573 6.158 0.000
conditionleading 0.087 0.014 68.677 6.018 0.000
conditiontrailing 0.028 0.014 65.692 2.013 0.048

Now add number of columns. Doesn’t converge with number of columns in random effects. This is also strong and clear, and the interactions suggest that all of this gets essentially canceled out in the three-column condition.

kable(summary(lmer(log(RT) ~ trial_num + condition * factor(n_col) + 
                     (condition | subnum), 
                   data = filter(dual_experts, 
                                 search_correct == 1, 
                                 abacus_correct == 1)))$coefficients, digits = 3)
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 1.201 0.035 75.790 34.014 0.000
trial_num -0.003 0.000 7337.382 -22.346 0.000
conditionout of play 0.120 0.019 161.625 6.378 0.000
conditionleading 0.154 0.018 179.412 8.408 0.000
conditiontrailing 0.064 0.018 187.765 3.569 0.000
factor(n_col)3 0.253 0.017 7249.594 15.336 0.000
conditionout of play:factor(n_col)3 -0.059 0.023 7243.509 -2.536 0.011
conditionleading:factor(n_col)3 -0.135 0.023 7246.326 -5.775 0.000
conditiontrailing:factor(n_col)3 -0.073 0.023 7250.239 -3.106 0.002

Is there an order effect?

dual_experts$second_half <- dual_experts$trial_num > 64
kable(summary(lmer(log(RT) ~ trial_num + condition * factor(n_col) * second_half + 
                     (condition | subnum), 
                   data = filter(dual_experts, 
                                 search_correct == 1, 
                                 abacus_correct == 1)))$coefficients, digits = 3)
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 1.208 0.037 93.170 32.477 0.000
trial_num -0.003 0.000 7333.356 -12.595 0.000
conditionout of play 0.102 0.025 488.360 4.060 0.000
conditionleading 0.151 0.025 554.449 6.104 0.000
conditiontrailing 0.059 0.025 619.098 2.406 0.016
factor(n_col)3 0.254 0.023 7274.846 10.983 0.000
second_halfTRUE 0.029 0.028 7276.300 1.065 0.287
conditionout of play:factor(n_col)3 -0.039 0.033 7285.226 -1.172 0.241
conditionleading:factor(n_col)3 -0.123 0.033 7290.802 -3.729 0.000
conditiontrailing:factor(n_col)3 -0.050 0.033 7314.139 -1.513 0.130
conditionout of play:second_halfTRUE 0.037 0.033 7289.411 1.112 0.266
conditionleading:second_halfTRUE 0.007 0.033 7304.564 0.198 0.843
conditiontrailing:second_halfTRUE 0.009 0.033 7313.240 0.261 0.794
factor(n_col)3:second_halfTRUE -0.003 0.033 7148.447 -0.092 0.927
conditionout of play:factor(n_col)3:second_halfTRUE -0.040 0.047 7275.376 -0.858 0.391
conditionleading:factor(n_col)3:second_halfTRUE -0.022 0.047 7299.465 -0.470 0.638
conditiontrailing:factor(n_col)3:second_halfTRUE -0.044 0.047 7275.484 -0.930 0.353

2 Experiment 2: Experts with no dual task

2.1 Data prep

Read data.

single_experts_raw <- read.csv("data/Upright Single Task Expert Data w. subject IDs and order.csv")
names(single_experts_raw) <- c("bead_type", "condition", "number_requested", 
                             "X_pos","Y_pos","search_correct","RT",
                             "abacus_val","abacus_correct","n_col",
                             "subnum","outlier","trial_type",
                             "subid","n_conditions","order","subid_xexpts")
single_experts_raw %<>% 
  mutate(bead_type = factor(bead_type, 
                            levels = c(0,1), 
                            labels = c("heaven","earth")),
         condition = factor(condition, 
                            levels = c(1,2,3,4), 
                            labels = c("in play", "out of play", 
                                       "leading","trailing")))

Exclusions. Filter pilot participants.

pilot_subs <- single_experts_raw %>%
  group_by(subnum) %>%
  summarise(pilot = any(n_col == 0)) %>%
  filter(pilot) 
  
single_experts <- filter(single_experts_raw, 
                       !subnum %in% pilot_subs$subnum)
single_experts %<>% 
  group_by(subnum) %>%
  mutate(trial_num = 1:n())

Check to make sure we have a consistent number of trials, no training trials.

qplot(subnum, trial_num, data = single_experts)

All participants have full data.

RT exclusions. Note that there are a few 0 RTs. What’s the deal with these?

sum(single_experts$RT == 0)
## [1] 43
single_experts %<>% filter(RT > 0, 
                         !is.na(RT))

Again clip in log space.

qplot(log(RT), data = single_experts, 
      fill = log(RT) > mean(log(RT)) + 3*sd(log(RT)) |
        log(RT) < mean(log(RT)) - 3*sd(log(RT)))

Clip these.

lmean <- mean(log(single_experts$RT))
lsd <- sd(log(single_experts$RT))
single_experts$RT[log(single_experts$RT) > lmean + 3*lsd |
                  log(single_experts$RT) < lmean - 3*lsd] <- NA

Replot in linear space just to check.

qplot(RT, data = single_experts)

Looks good.

2.2 Demographics and descriptives

kable(single_experts %>% 
        group_by(subnum) %>% 
        summarise(n = n()) %>% 
        summarise(n = n()))

2.3 n

67

kable(single_experts %>% 
        group_by(subnum) %>% 
        summarise(RT = mean(RT, na.rm=TRUE)) %>% 
        mutate(condition = "all") %>% group_by(condition) %>%
        multi_boot_standard(col = "RT"))
condition mean ci_lower ci_upper
all 1.969064 1.847767 2.100203

2.4 RT and accuracy analyses

Basic analyses.

ms <- single_experts %>%
  filter(search_correct == 1) %>%
  group_by(subnum, condition) %>%
  summarise(RT = mean(RT, na.rm=TRUE)) %>%
  group_by(condition) %>%
  multi_boot_standard(col = "RT")

kable(ms, digits = 2)
condition mean ci_lower ci_upper
in play 1.88 1.76 2.01
out of play 2.09 1.96 2.22
leading 2.05 1.91 2.18
trailing 1.89 1.76 2.01 
ggplot(ms,aes(x = condition, y = mean, fill = condition)) + 
  geom_bar(stat = "identity") + 
  geom_linerange(aes(ymin = ci_lower, ymax = ci_upper))

Add two other variables: bead type and number of columns.

ms <- single_experts %>%
  filter(search_correct == 1) %>%
  group_by(subnum, condition, bead_type, n_col) %>%
  summarise(RT = mean(RT, na.rm=TRUE)) %>%
  group_by(condition, bead_type, n_col) %>%
  multi_boot_standard(col = "RT")

kable(ms, digits = 2)
condition bead_type n_col mean ci_lower ci_upper
in play heaven 2 1.78 1.66 1.92
in play heaven 3 2.10 1.92 2.28
in play earth 2 1.68 1.57 1.79
in play earth 3 2.04 1.88 2.20
out of play heaven 2 2.14 2.00 2.29
out of play heaven 3 2.48 2.30 2.67
out of play earth 2 1.82 1.70 1.94
out of play earth 3 2.22 2.05 2.39
leading heaven 2 2.10 1.94 2.29
leading heaven 3 2.57 2.37 2.78
leading earth 2 1.80 1.68 1.93
leading earth 3 2.09 1.93 2.26
trailing heaven 2 2.00 1.85 2.18
trailing heaven 3 2.39 2.20 2.59
trailing earth 2 1.64 1.53 1.75
trailing earth 3 1.94 1.80 2.10 
ggplot(ms, aes(x = condition, y = mean, fill = condition)) + 
  geom_bar(stat = "identity") + 
  geom_linerange(aes(ymin = ci_lower, ymax = ci_upper)) + 
  facet_grid(bead_type ~ n_col)

Now in this experiment, we’re actually seeing effects in the three-column case. That suggests that it was the load of the three-column abacus reading that was suppressing the attentional effects, which is actually kind of nice and interesting.

2.5 Stats

Same LMER as before. This time the model didn’t converge with random condition effects.

kable(summary(lmer(log(RT) ~ trial_num + condition + 
                     (1 | subnum), 
                   data = filter(single_experts, 
                                 search_correct == 1)))$coefficients, 
      digits = 3)
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 0.533 0.033 86.563 16.301 0.000
trial_num -0.001 0.000 8045.513 -4.481 0.000
conditionout of play 0.118 0.014 8044.111 8.729 0.000
conditionleading 0.085 0.014 8044.108 6.276 0.000
conditiontrailing 0.013 0.013 8044.112 0.928 0.353

Now add number of columns again. Here there are no interactions, which is clear and nice. Interestingly, now the model will converge with condition in the random effects.

kable(summary(lmer(log(RT) ~ trial_num + condition * factor(n_col) + 
                     (condition | subnum), 
                   data = filter(single_experts, 
                                 search_correct == 1)))$coefficients, digits = 3)
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 0.452 0.035 85.814 12.792 0.000
trial_num -0.001 0.000 8028.483 -4.533 0.000
conditionout of play 0.119 0.020 209.993 6.011 0.000
conditionleading 0.092 0.019 401.992 4.819 0.000
conditiontrailing 0.019 0.020 202.561 0.953 0.342
factor(n_col)3 0.163 0.019 7915.599 8.705 0.000
conditionout of play:factor(n_col)3 -0.003 0.027 7915.755 -0.107 0.914
conditionleading:factor(n_col)3 -0.016 0.027 7914.997 -0.599 0.549
conditiontrailing:factor(n_col)3 -0.014 0.026 7918.387 -0.535 0.593

3 Experiments 1 and 2 together

Bind everything together.

experts <- bind_rows(filter(single_experts,
                            search_correct ==1) %>%
                       mutate(expt = "single task", 
                              group = "experts"),
                     filter(dual_experts, 
                            search_correct == 1, 
                            abacus_correct == 1) %>%
                       mutate(expt = "dual task", 
                              group = "experts")) %>%
    mutate(subid_xexpts = as.numeric(as.character(subid_xexpts)), 
           order = ifelse(order == "1", 1, ifelse(order == "N/A", NA, 2)))

3.1 Visualization

We don’t learn much more this way, but we can plot everything together

ms <- experts %>%
  group_by(subnum, condition, expt) %>%
  summarise(RT = mean(RT, na.rm=TRUE)) %>%
  group_by(condition, expt) %>%
  multi_boot_standard(col = "RT")

ggplot(ms,aes(x = condition, y = mean, fill = condition)) + 
  geom_bar(stat = "identity") + 
  geom_linerange(aes(ymin = ci_lower, ymax = ci_upper)) + 
  facet_grid(~expt)

3.2 Stats

Model the whole thing together.

kable(summary(lmer(log(RT) ~ expt * trial_num + expt * condition + 
                     (condition | subnum), 
                   data = filter(experts)))$coefficients, 
      digits = 3)
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 1.330 0.035 142.996 38.151 0.000
exptsingle task -0.797 0.049 142.599 -16.416 0.000
trial_num -0.003 0.000 15381.100 -20.279 0.000
conditionout of play 0.092 0.015 133.803 6.064 0.000
conditionleading 0.087 0.014 156.398 6.080 0.000
conditiontrailing 0.028 0.015 138.986 1.887 0.061
exptsingle task:trial_num 0.002 0.000 15383.058 11.321 0.000
exptsingle task:conditionout of play 0.026 0.021 131.410 1.259 0.210
exptsingle task:conditionleading -0.002 0.020 154.793 -0.087 0.931
exptsingle task:conditiontrailing -0.016 0.021 135.491 -0.760 0.449

This is very interpretable:

  • the single task is faster than dual,
  • you get faster as you continue doing the tasks
  • you don’t improve as much with repeated trials in the single task (because you’re not getting faster at reading the abacus)
  • there are effects of in vs. out of play and leading vs. following
  • there are no interactions with task.

Awesome.

3.3 Order effects

First figure out who participated in both. You only got a subid_xexpt if you participated in both?

xexpts <- unique(experts$subid_xexpts)
length(xexpts[!is.na(xexpts)])
## [1] 17

Next, check the stats for this group. First same model as above.

kable(summary(lmer(log(RT) ~ expt * trial_num + expt * condition + 
                     (condition | subnum), 
                   data = filter(experts, !is.na(subid_xexpts))))$coefficients, 
      digits = 3)
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 1.160 0.073 34.815 15.989 0.000
exptsingle task -0.663 0.103 34.843 -6.459 0.000
trial_num -0.003 0.000 4032.570 -10.684 0.000
conditionout of play 0.122 0.029 32.523 4.209 0.000
conditionleading 0.115 0.027 52.434 4.266 0.000
conditiontrailing 0.061 0.028 40.645 2.152 0.037
exptsingle task:trial_num 0.002 0.000 4033.188 5.102 0.000
exptsingle task:conditionout of play 0.008 0.041 32.214 0.208 0.837
exptsingle task:conditionleading -0.029 0.038 52.243 -0.746 0.459
exptsingle task:conditiontrailing -0.048 0.040 39.763 -1.220 0.230

Next, add interactions of order.

print("foo")
## [1] "foo"
kable(summary(lmer(log(RT) ~ expt * trial_num + expt * condition * order + 
                     (condition | subnum), 
                   data = filter(experts, !is.na(subid_xexpts))))$coefficients, 
      digits = 3)
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 1.306 0.282 30.338 4.633 0.000
exptsingle task -0.675 0.357 30.317 -1.894 0.068
trial_num -0.003 0.000 4027.067 -10.681 0.000
conditionout of play 0.173 0.117 32.449 1.477 0.149
conditionleading 0.121 0.109 45.933 1.111 0.272
conditiontrailing 0.186 0.111 40.012 1.685 0.100
order -0.086 0.159 30.156 -0.537 0.595
exptsingle task:trial_num 0.002 0.000 4027.813 5.090 0.000
exptsingle task:conditionout of play -0.081 0.147 31.484 -0.555 0.583
exptsingle task:conditionleading 0.012 0.137 44.967 0.087 0.931
exptsingle task:conditiontrailing -0.116 0.139 38.774 -0.837 0.408
exptsingle task:order -0.013 0.220 30.116 -0.057 0.955
conditionout of play:order -0.030 0.066 31.854 -0.450 0.656
conditionleading:order -0.003 0.062 45.218 -0.053 0.958
conditiontrailing:order -0.073 0.062 39.456 -1.176 0.247
exptsingle task:conditionout of play:order 0.058 0.090 31.166 0.643 0.525
exptsingle task:conditionleading:order -0.031 0.085 44.773 -0.370 0.713
exptsingle task:conditiontrailing:order 0.030 0.086 38.591 0.353 0.726

Hard to conclude anything from this. Let’s look at the plot.

ms <- experts %>%
  filter(!is.na(subid_xexpts)) %>%
  group_by(subnum, condition, expt, order) %>%
  summarise(RT = mean(RT, na.rm=TRUE)) %>%
  group_by(condition, expt, order) %>%
  multi_boot_standard(col = "RT")

ggplot(ms,aes(x = condition, y = mean, fill = condition)) + 
  geom_bar(stat = "identity") + 
  geom_linerange(aes(ymin = ci_lower, ymax = ci_upper)) + 
  facet_grid(order~expt)

4 Experiment 3: Naive participants, single task

4.1 Data prep

Read data.

naive_raw <- read.csv("data/UprightAdultData.csv")
names(naive_raw) <- c("bead_type", "condition", "number_requested", 
                             "X_pos","Y_pos","search_correct","RT",
                             "abacus_val","abacus_correct","n_col",
                             "subnum","outlier","trial_type", "expertise")
naive_raw %<>% 
  mutate(bead_type = factor(bead_type, 
                            levels = c(0,1), 
                            labels = c("heaven","earth")),
         condition = factor(condition, 
                            levels = c(1,2,3,4), 
                            labels = c("in play", "out of play", 
                                       "leading","trailing")))

Exclusions. Filter pilot participants.

pilot_subs <- naive_raw %>%
  group_by(subnum) %>%
  summarise(pilot = any(n_col == 0)) %>%
  filter(pilot) 
  
naive <- filter(naive_raw, 
                !subnum %in% pilot_subs$subnum)
naive %<>% 
  group_by(subnum) %>%
  mutate(trial_num = 1:n())

Check to make sure we have a consistent number of trials, no training trials.

qplot(subnum, trial_num, data = naive)

All participants have full data.

RT exclusions. Again clip in log space.

qplot(log(RT), data = naive, 
      fill = log(RT) > mean(log(RT)) + 3*sd(log(RT)) |
        log(RT) < mean(log(RT)) - 3*sd(log(RT)))

Clip these.

lmean <- mean(log(naive$RT))
lsd <- sd(log(naive$RT))
naive$RT[log(naive$RT) > lmean + 3*lsd |
           log(naive$RT) < lmean - 3*lsd] <- NA

Replot in linear space just to check.

qplot(RT, data = naive)

Looks good.

4.2 Demographics and descriptives

kable(naive %>% 
        group_by(subnum) %>% 
        summarise(n = n()) %>% 
        summarise(n = n()))

4.3 n

56

kable(naive %>% 
        group_by(subnum) %>% 
        summarise(RT = mean(RT, na.rm=TRUE)) %>% 
        mutate(condition = "all") %>% group_by(condition) %>%
        multi_boot_standard(col = "RT"))
condition mean ci_lower ci_upper
all 1.195807 1.11825 1.270167

4.4 RT and accuracy analyses

Basic analyses. Summary: In this experiment, the effects are smaller, but still present, at all column levels.

ms <- naive %>%
  filter(search_correct == 1) %>%
  group_by(subnum, condition) %>%
  summarise(RT = mean(RT, na.rm=TRUE)) %>%
  group_by(condition) %>%
  multi_boot_standard(col = "RT")

kable(ms, digits = 2)
condition mean ci_lower ci_upper
in play 1.15 1.07 1.24
out of play 1.23 1.15 1.30
leading 1.23 1.15 1.31
trailing 1.16 1.10 1.24 
ggplot(ms,aes(x = condition, y = mean, fill = condition)) + 
  geom_bar(stat = "identity") + 
  geom_linerange(aes(ymin = ci_lower, ymax = ci_upper))

Add two other variables: bead type and number of columns.

ms <- naive %>%
  filter(search_correct == 1) %>%
  group_by(subnum, condition, bead_type, n_col) %>%
  summarise(RT = mean(RT, na.rm=TRUE)) %>%
  group_by(condition, bead_type, n_col) %>%
  multi_boot_standard(col = "RT")

kable(ms, digits = 2)
condition bead_type n_col mean ci_lower ci_upper
in play heaven 2 1.09 1.01 1.17
in play heaven 3 1.22 1.11 1.33
in play earth 2 1.10 1.02 1.20
in play earth 3 1.20 1.11 1.30
out of play heaven 2 1.22 1.13 1.31
out of play heaven 3 1.35 1.24 1.46
out of play earth 2 1.14 1.07 1.21
out of play earth 3 1.29 1.21 1.38
leading heaven 2 1.22 1.13 1.31
leading heaven 3 1.35 1.24 1.47
leading earth 2 1.14 1.07 1.22
leading earth 3 1.29 1.20 1.38
trailing heaven 2 1.19 1.11 1.29
trailing heaven 3 1.33 1.22 1.45
trailing earth 2 1.08 1.00 1.17
trailing earth 3 1.19 1.12 1.27 
ggplot(ms, aes(x = condition, y = mean, fill = condition)) + 
  geom_bar(stat = "identity") + 
  geom_linerange(aes(ymin = ci_lower, ymax = ci_upper)) + 
  facet_grid(bead_type ~ n_col)

Add abacus expertise instead.

ms <- naive %>%
  filter(search_correct == 1) %>%
  group_by(subnum, condition, expertise) %>%
  summarise(RT = mean(RT, na.rm=TRUE)) %>%
  group_by(condition, expertise) %>%
  multi_boot_standard(col = "RT")

kable(ms, digits = 2)
condition expertise mean ci_lower ci_upper
in play 0 1.20 1.08 1.32
in play 1 1.06 0.96 1.17
in play 2 1.21 0.92 1.55
out of play 0 1.27 1.16 1.38
out of play 1 1.14 1.05 1.25
out of play 2 1.30 1.12 1.48
leading 0 1.25 1.14 1.36
leading 1 1.14 1.03 1.26
leading 2 1.40 1.18 1.64
trailing 0 1.21 1.11 1.32
trailing 1 1.07 0.99 1.16
trailing 2 1.22 1.01 1.53 
ggplot(ms, aes(x = condition, y = mean, fill = condition)) + 
  geom_bar(stat = "identity") + 
  geom_linerange(aes(ymin = ci_lower, ymax = ci_upper)) + 
  facet_grid(. ~ expertise)

4.5 Stats

Same LMER as in the previous two, with the full random effects structure. This model shows both effects, as before.

kable(summary(lmer(log(RT) ~ trial_num + condition + 
                     (condition | subnum), 
                   data = filter(naive, 
                                 search_correct == 1)))$coefficients, 
      digits = 3)
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 0.118 0.034 59.745 3.459 0.001
trial_num -0.001 0.000 6891.010 -10.368 0.000
conditionout of play 0.081 0.012 133.951 6.520 0.000
conditionleading 0.077 0.014 64.863 5.671 0.000
conditiontrailing 0.023 0.013 107.088 1.764 0.081

Now add number of columns again. Here again there are no interactions, suggesting that column didn’t affect matters.

kable(summary(lmer(log(RT) ~ trial_num + condition * factor(n_col) + 
                     (condition | subnum), 
                   data = filter(naive, 
                                 search_correct == 1)))$coefficients, digits = 3)
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 0.083 0.035 66.218 2.367 0.021
trial_num -0.001 0.000 6885.169 -10.529 0.000
conditionout of play 0.067 0.017 405.167 4.048 0.000
conditionleading 0.063 0.018 153.058 3.567 0.000
conditiontrailing 0.011 0.017 181.253 0.676 0.500
factor(n_col)3 0.071 0.016 6794.703 4.542 0.000
conditionout of play:factor(n_col)3 0.029 0.022 6792.594 1.292 0.196
conditionleading:factor(n_col)3 0.030 0.022 6795.768 1.350 0.177
conditiontrailing:factor(n_col)3 0.022 0.022 6795.111 1.011 0.312

And check for interactions of expertise level. First with expertise as a continuous variable.

kable(summary(lmer(log(RT) ~ trial_num + condition * expertise + 
                     (condition | subnum), 
                   data = filter(naive, 
                                 search_correct == 1)))$coefficients, digits = 3)
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 0.136 0.043 56.867 3.138 0.003
trial_num -0.001 0.000 6889.146 -10.387 0.000
conditionout of play 0.075 0.016 134.303 4.712 0.000
conditionleading 0.051 0.017 70.439 3.093 0.003
conditiontrailing 0.022 0.016 95.912 1.339 0.184
expertise -0.033 0.049 54.030 -0.670 0.505
conditionout of play:expertise 0.013 0.018 134.600 0.686 0.494
conditionleading:expertise 0.049 0.019 70.380 2.535 0.013
conditiontrailing:expertise 0.001 0.019 95.103 0.054 0.957

Now with expertise as a factor.

kable(summary(lmer(log(RT) ~ trial_num + condition * factor(expertise) + 
                     (condition | subnum), 
                   data = filter(naive, 
                                 search_correct == 1)))$coefficients, digits = 3)
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 0.152 0.045 55.637 3.398 0.001
trial_num -0.001 0.000 6886.050 -10.400 0.000
conditionout of play 0.077 0.017 126.736 4.655 0.000
conditionleading 0.056 0.017 70.759 3.206 0.002
conditiontrailing 0.023 0.017 88.689 1.323 0.189
factor(expertise)1 -0.105 0.073 53.019 -1.435 0.157
factor(expertise)2 0.003 0.111 53.036 0.024 0.981
conditionout of play:factor(expertise)1 0.000 0.028 125.384 -0.005 0.996
conditionleading:factor(expertise)1 0.029 0.029 70.489 1.011 0.315
conditiontrailing:factor(expertise)1 -0.003 0.029 87.583 -0.107 0.915
conditionout of play:factor(expertise)2 0.037 0.042 127.671 0.895 0.372
conditionleading:factor(expertise)2 0.116 0.044 70.735 2.651 0.010
conditiontrailing:factor(expertise)2 0.006 0.043 87.705 0.140 0.889

In both cases, we see a slightly bigger effect of leading zeros being slower for the participants with more abacus exposure. This is not totally unreasonable, but it’s a small effect.

Finally, check a followup model for only zero expertise subjects. Does not converge with random slope.

kable(summary(lmer(log(RT) ~ trial_num + condition + 
                     (1 | subnum), 
                   data = filter(naive, 
                                 expertise == 0,
                                 search_correct == 1)))$coefficients, digits = 3)
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 0.164 0.043 37.349 3.833 0.000
trial_num -0.001 0.000 3930.019 -9.279 0.000
conditionout of play 0.078 0.015 3930.005 5.264 0.000
conditionleading 0.056 0.015 3930.008 3.768 0.000
conditiontrailing 0.023 0.015 3930.013 1.570 0.116

5 Experiments 2 and 3

Bind all data.

d <- bind_rows(filter(naive, 
                      search_correct == 1) %>%
                 mutate(expt = "single task",
                        group = "naive"), 
               experts)

Now plot.

ms <- d %>%
  filter(expt == "single task") %>%
  group_by(subnum, condition, group) %>%
  summarise(RT = mean(RT, na.rm=TRUE)) %>%
  group_by(condition, group) %>%
  multi_boot_standard(col = "RT")

ggplot(ms, aes(x = condition, y = mean, fill = condition)) + 
  geom_bar(stat = "identity") + 
  geom_linerange(aes(ymin = ci_lower, ymax = ci_upper)) + 
  facet_grid( ~ group)

And a statistical model.

kable(summary(lmer(log(RT) ~ trial_num + condition * group + 
                     (condition | subnum), 
                   data = filter(d, 
                                 expt == "single task")))$coefficients, 
      digits = 3)
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 0.549 0.032 128.397 16.937 0.000
trial_num -0.001 0.000 14930.532 -9.702 0.000
conditionout of play 0.118 0.013 163.223 9.056 0.000
conditionleading 0.085 0.013 133.027 6.493 0.000
conditiontrailing 0.013 0.014 124.919 0.930 0.354
groupnaive -0.450 0.047 120.917 -9.499 0.000
conditionout of play:groupnaive -0.037 0.019 158.820 -1.950 0.053
conditionleading:groupnaive -0.009 0.019 129.612 -0.444 0.658
conditiontrailing:groupnaive 0.009 0.020 121.782 0.469 0.640

Summary - experts are way slower than naive participants, even when they don’t have the dual task. That’s interesting, I think - they are still processing the abacus somehow prior to searching.

There’s not much in the way of interactions of the effects with group, with one exception. The naive participants don’t show as big an “out of play” effect as the experts, trending with p = 0.0511761. So there may be some small difference there.

6 All data

6.1 Plots for paper

ms <- d %>%
  group_by(subnum, condition, expt, expertise, group) %>%
  summarise(RT = mean(RT, na.rm=TRUE)) %>%
  group_by(condition, expt, expertise, group) %>%
  multi_boot_standard(col = "RT")

First in-play/out-of-play.

ggplot(filter(ms, group=="experts", 
              condition == "in play" | condition == "out of play"), 
       aes(x = condition, y = mean, fill = condition)) + 
  geom_bar(stat = "identity") + 
  geom_linerange(aes(ymin = ci_lower, ymax = ci_upper)) + 
  facet_grid( ~ expt) + 
  theme_mikabr() + 
  scale_fill_solarized()

Then trailing/leading.

ggplot(filter(ms, group=="experts", 
              condition == "leading" | condition == "trailing"), 
       aes(x = condition, y = mean, fill = condition)) + 
  geom_bar(stat = "identity") + 
  geom_linerange(aes(ymin = ci_lower, ymax = ci_upper)) + 
  facet_grid( ~ expt) + 
  theme_mikabr() + 
  scale_fill_solarized()

Then the same in/out of play by expertise.

ggplot(filter(ms, group=="naive", 
              condition == "in play" | condition == "out of play"), 
       aes(x = condition, y = mean, fill = condition)) + 
  geom_bar(stat = "identity") + 
  geom_linerange(aes(ymin = ci_lower, ymax = ci_upper)) + 
  facet_grid( ~ expertise) + 
  theme_mikabr() + 
  scale_fill_solarized()

and leading/trailing by expertise.

ms$expertise <- factor(ms$expertise, levels = c(0, 1, 2), labels = c("No experience","Recognition only","Experience using"))
ggplot(filter(ms, group=="naive", 
              condition == "leading" | condition == "trailing"), 
       aes(x = condition, y = mean, fill = condition)) + 
  geom_bar(stat = "identity") + 
  geom_linerange(aes(ymin = ci_lower, ymax = ci_upper)) + 
  facet_grid( ~ expertise) + 
  theme_mikabr() + 
  scale_fill_solarized()

6.2 Basic stats

Consolidate all data into a single model.

kable(summary(lmer(log(RT) ~ trial_num * group * expt + 
                     condition * group * expt + 
                     (condition | subnum), 
                   data = d))$coefficients, 
      digits = 3)
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 1.330 0.034 203.446 39.045 0.000
trial_num -0.003 0.000 22271.939 -21.355 0.000
groupnaive -0.416 0.049 201.474 -8.514 0.000
exptsingle task -0.797 0.047 202.911 -16.801 0.000
conditionout of play 0.092 0.014 214.582 6.555 0.000
conditionleading 0.087 0.014 207.265 6.248 0.000
conditiontrailing 0.028 0.014 198.718 1.992 0.048
trial_num:groupnaive -0.001 0.000 22278.563 -3.150 0.002
trial_num:exptsingle task 0.002 0.000 22273.093 11.929 0.000
groupnaive:conditionout of play -0.037 0.020 201.129 -1.860 0.064
groupnaive:conditionleading -0.008 0.020 198.235 -0.391 0.696
groupnaive:conditiontrailing 0.010 0.020 183.728 0.499 0.618
exptsingle task:conditionout of play 0.026 0.019 210.530 1.364 0.174
exptsingle task:conditionleading -0.002 0.019 205.236 -0.089 0.929
exptsingle task:conditiontrailing -0.016 0.020 193.654 -0.805 0.422

In the full model we see all the effects holding up, with not much evidence for interactions of expertise or experiment type. That’s pretty much exactly what we thought was going on.

6.3 Post-hoc analysis: Targets closer to the beam

Does position mediate the “in play/out of play” effect?

First, the baseline model for this analysis: only earthly beads (since heavenly bead position and “in play” is confounded.

There’s an interesting decision to make here in the random effects. I think it’s best to go with no condition or Y_pos random effect both for convergence reasons and ebcause it’s hard theoretically to interpret when you have all the interactions here.

kable(summary(lmer(log(RT) ~ trial_num * group * expt + 
                     condition * group * expt + 
                     (1 | subnum), 
                   data = filter(d, 
                                 condition %in% c("in play","out of play"),
                                 bead_type == "earth")))$coefficients, 
      digits = 3)
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 1.300 0.035 270.637 37.549 0.000
trial_num -0.002 0.000 8269.966 -12.262 0.000
groupnaive -0.400 0.049 264.586 -8.083 0.000
exptsingle task -0.784 0.048 268.799 -16.288 0.000
conditionout of play 0.091 0.014 8255.969 6.263 0.000
trial_num:groupnaive -0.001 0.000 8264.927 -2.063 0.039
trial_num:exptsingle task 0.002 0.000 8268.280 7.002 0.000
groupnaive:conditionout of play -0.020 0.020 8254.536 -0.971 0.331
exptsingle task:conditionout of play -0.001 0.020 8255.449 -0.025 0.980

Next, add vertical position as a predictor.

kable(summary(lmer(log(RT) ~ trial_num * group * expt + 
                     condition * group * expt + 
                     Y_pos * group * expt + 
                     (1 | subnum), 
                   data = filter(d, 
                                 condition %in% c("in play","out of play"),
                                 bead_type == "earth")))$coefficients, 
      digits = 3)
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 1.177 0.045 757.722 26.232 0.000
trial_num -0.002 0.000 8267.064 -12.377 0.000
groupnaive -0.457 0.064 723.261 -7.178 0.000
exptsingle task -0.790 0.062 743.118 -12.712 0.000
conditionout of play 0.019 0.022 8266.453 0.842 0.400
Y_pos 0.032 0.007 8276.033 4.294 0.000
trial_num:groupnaive -0.001 0.000 8261.886 -2.083 0.037
trial_num:exptsingle task 0.002 0.000 8265.250 6.958 0.000
groupnaive:conditionout of play -0.056 0.032 8260.536 -1.764 0.078
exptsingle task:conditionout of play -0.006 0.031 8262.576 -0.194 0.846
groupnaive:Y_pos 0.015 0.011 8266.453 1.412 0.158
exptsingle task:Y_pos 0.002 0.010 8269.996 0.196 0.845

Now we see a significant effect of Y position with no remaining effect of “out of play.” Check for interactions?

kable(summary(lmer(log(RT) ~ trial_num * group * expt + 
                     Y_pos * group * expt * condition + 
                     (1 | subnum), 
                   data = filter(d, 
                                 condition %in% c("in play","out of play"),
                                 bead_type == "earth")))$coefficients, 
      digits = 3)
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 1.186 0.053 1451.318 22.177 0.000
trial_num -0.002 0.000 8264.037 -12.389 0.000
groupnaive -0.554 0.076 1395.006 -7.292 0.000
exptsingle task -0.684 0.074 1435.486 -9.215 0.000
Y_pos 0.030 0.011 8272.173 2.795 0.005
conditionout of play -0.005 0.077 8270.315 -0.070 0.944
trial_num:groupnaive -0.001 0.000 8258.967 -2.093 0.036
trial_num:exptsingle task 0.002 0.000 8262.229 6.953 0.000
groupnaive:Y_pos 0.040 0.015 8264.475 2.671 0.008
exptsingle task:Y_pos -0.025 0.015 8268.709 -1.719 0.086
Y_pos:conditionout of play 0.005 0.015 8270.444 0.324 0.746
groupnaive:conditionout of play 0.188 0.111 8264.227 1.696 0.090
exptsingle task:conditionout of play -0.273 0.108 8268.957 -2.535 0.011
groupnaive:Y_pos:conditionout of play -0.049 0.021 8264.941 -2.310 0.021
exptsingle task:Y_pos:conditionout of play 0.054 0.021 8269.564 2.591 0.010

Interestingly, we are seeing some three-way interactions, but these are a bit hard to interpret. Let’s subset to Experiments 1, 2, and 3.

kable(summary(lmer(log(RT) ~ trial_num + 
                     Y_pos * condition + 
                     (1 | subnum), 
                   data = filter(d, 
                                 group == "experts", 
                                 expt == "dual task",
                                 condition %in% c("in play","out of play"),
                                 bead_type == "earth")))$coefficients, 
      digits = 3)
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 1.186 0.053 493.180 22.201 0.000
trial_num -0.002 0.000 2727.389 -12.390 0.000
Y_pos 0.030 0.011 2730.094 2.796 0.005
conditionout of play -0.005 0.077 2729.476 -0.070 0.944
Y_pos:conditionout of play 0.005 0.015 2729.519 0.324 0.746 
kable(summary(lmer(log(RT) ~ trial_num + 
                     Y_pos * condition + 
                     (1 | subnum), 
                   data = filter(d, 
                                 group == "experts", 
                                 expt == "single task",
                                 condition %in% c("in play","out of play"),
                                 bead_type == "earth")))$coefficients, 
      digits = 3)
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 0.502 0.057 580.806 8.882 0.000
trial_num -0.001 0.000 2967.937 -2.613 0.009
Y_pos 0.004 0.011 2969.837 0.366 0.714
conditionout of play -0.279 0.084 2970.836 -3.336 0.001
Y_pos:conditionout of play 0.059 0.016 2971.285 3.639 0.000 
kable(summary(lmer(log(RT) ~ trial_num + 
                     Y_pos * condition + 
                     (1 | subnum), 
                   data = filter(d, 
                                 group == "naive", 
                                 expt == "single task",
                                 condition %in% c("in play","out of play"),
                                 bead_type == "earth")))$coefficients, 
      digits = 3)
Estimate Std. Error df t value Pr(>|t|)
(Intercept) -0.052 0.049 333.905 -1.056 0.292
trial_num -0.001 0.000 2565.454 -6.595 0.000
Y_pos 0.044 0.009 2567.042 4.763 0.000
conditionout of play -0.090 0.069 2566.383 -1.308 0.191
Y_pos:conditionout of play 0.009 0.013 2566.463 0.715 0.475

So we see the same unpredicted effect in Experiment 2 as in the previous analysis.