Script B: Experiment 1

In [1]:

library(osfr)
library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.1     ✔ tibble    3.2.1
✔ lubridate 1.9.3     ✔ tidyr     1.3.1
✔ purrr     1.0.2     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(knitr)
library(gt)
library(gtExtras)
library(rstatix)

Attaching package: 'rstatix'

The following object is masked from 'package:stats':

    filter

library(ggpubr)
library(here)

here() starts at /Users/rp3650/Library/CloudStorage/GoogleDrive-robpetrosino@gmail.com/My Drive/Academics/projects/morphology/morphological-decomposition/sub-projects/frequency-effects/frequency-effect_masked-priming

library(lme4)

Loading required package: Matrix

Attaching package: 'Matrix'

The following objects are masked from 'package:tidyr':

    expand, pack, unpack

library(car)

Loading required package: carData

Attaching package: 'car'

The following object is masked from 'package:dplyr':

    recode

The following object is masked from 'package:purrr':

    some

library(emmeans)

Welcome to emmeans.
Caution: You lose important information if you filter this package's results.
See '? untidy'

Materials

In [2]:

#|
words <- read_csv('../materials/experiments-1-2/frequency-effects_experiments-1-2_stimuli.csv') %>%
            filter(type=='word')

Rows: 208 Columns: 15
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr  (3): Word, type, freq.bin
dbl (12): Mean_RT, Mean_Accuracy, SUBTLWF, Freq_HAL, Length, Ortho_N, BG_Sum...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

In [3]:

#|
needs_mutated <- function(x) {
  if(!is.numeric(x)) return(FALSE)
  max(x, na.rm = TRUE) > 100
}

words %>%
  group_by(freq.bin) %>%
  mutate(Freq_HAL.Pm = (Freq_HAL * 10^6)/(131 * 10^6)) %>%
  summarise(N = n(),
            minFreq= min(Freq_HAL.Pm, na.rm = T), maxFreq=max(Freq_HAL.Pm, na.rm = T),
            meanFreq = mean(Freq_HAL.Pm, na.rm = T), sdFreq = sd(Freq_HAL.Pm, na.rm = T),
            minSUBFreq= min(SUBTLWF, na.rm = T), maxSUBFreq=max(SUBTLWF, na.rm = T),
            meanSUBFreq = mean(SUBTLWF, na.rm = T), sdSUBFreq = sd(SUBTLWF, na.rm = T),
            #meanRT=mean(as.numeric(I_Mean_RT), na.rm = T), sdRT=sd(as.numeric(I_Mean_RT), na.rm = T),
            #meanLength = mean(Length, na.rm = T), sdLength=sd(Length, na.rm = T)
            ) %>% 
  mutate_if(is.numeric, ~round(.)) %>% 
  mutate_if(needs_mutated, ~round(.)) %>%
  gt() %>%
  cols_label(
    freq.bin = "frequency", 
    minFreq = "min", maxFreq = "max", meanFreq = "mean", sdFreq = "SD", 
    minSUBFreq = "min", maxSUBFreq = "max", meanSUBFreq = "mean", sdSUBFreq = "SD") %>%
  tab_spanner(
    label=md("**HAL**"),
    columns=3:6
  ) %>%
  tab_spanner(
    label=md("**SUBTLEX~US~**"),
    columns=7:10
  )

Descriptive statistics of the word items used. For both frequency databases, the word frequencies were converted to per-million count to ensure cross-comparison.

frequency	N	HAL				SUBTLEXUS
		min	max	mean	SD	min	max	mean	SD
high	52	45	4984	573	808	57	2691	210	388
low	52	6	570	64	93	7	24	13	5

Data analysis

In [4]:

# load the raw data dataframe
exp1_data_folder <- "data/experiment1"
exp1_rawdata_filename <- "experiment_1_preprocessed_data.csv"

## 02. check if the rawdata file exists. if not, download it from OSF.
if (!file.exists(here(exp1_data_folder, exp1_rawdata_filename))) {
  osf_retrieve_file("vn3r2") |> 
    osf_download(path = here(exp1_data_folder),
                 conflicts = "overwrite") 
}

## 03. read the data into R.
exp1_rawdata <- here(exp1_data_folder, exp1_rawdata_filename) |>
  read.csv(na = c("", "NA")) %>%
  mutate(primeTime = primeDuration - maskDuration) # calculating the actual SOA

exp1_info <- list()
exp1_info$intended_prime_duration <- 33
exp1_info$prime_dur_lb <- 25
exp1_info$prime_dur_ub <- 60
exp1_info$rt_lb <- 200
exp1_info$rt_ub <- 1800
exp1_info$freq_conditions <- c("high", "low", "non-word")
exp1_info$n_recruited <- exp1_rawdata$Rec_Session_Id |>
  unique() |>
  length()

exp1_rawdata.sub <- exp1_rawdata %>%
  filter(!is.na(TimeMeasure_Mean) & !is.na(primeDuration) & !is.na(responseError))
            
exp1_info$summary <- with(
  transform(exp1_rawdata.sub,
    RT_inrange = ifelse(RT >= exp1_info$rt_lb & RT <= exp1_info$rt_ub, 1, 0),
    Prime_inrange = ifelse((primeDuration - maskDuration) >= exp1_info$prime_dur_lb &
                             (primeDuration - maskDuration) <= exp1_info$prime_dur_ub, 1, 0)),
  {
    data.frame(aggregate(Start_Time ~ Rec_Session_Id + Crowdsourcing_SubjId, data=exp1_rawdata.sub, unique),
               aggregate(End_Time_Local ~ Rec_Session_Id + Crowdsourcing_SubjId, data=exp1_rawdata.sub, unique),
              aggregate(cbind(list, SelectedGender, SelectedAge, TimeMeasure_Mean, TimeMeasure_Std) ~ Rec_Session_Id + Crowdsourcing_SubjId, data=exp1_rawdata.sub, unique),
      aggregate(cbind(responseError, RT_inrange, Prime_inrange) ~ Rec_Session_Id + Crowdsourcing_SubjId, mean, data=exp1_rawdata.sub)
  )
}
)

exp1_info$summary <- exp1_info$summary[, -grep("Rec_Session_Id.|Crowdsourcing_SubjId.", colnames(exp1_info$summary))] # remove all extra aggregating columns (subj ID)

exp1_info$summary$Duration <- interval(ymd_hms(exp1_info$summary$Start_Time), 
                                             ymd_hms(exp1_info$summary$End_Time_Local)) |>
                                      lapply(function(interval_value) {interval_value/dminutes(1)}) |> 
                                           unlist()

Step 1: subject and item performance

In [5]:

exp1_step1_goodsubj <- exp1_info$summary |>
  subset(responseError <= .3) 

exp1_step1_subj_remain <- exp1_step1_goodsubj |> nrow()

exp1_step1_item.err <- exp1_rawdata.sub %>% group_by(condition_rec, target_rec) %>%
  summarise(word.percent=mean(responseError)*100) %>% 
  filter(word.percent > 30)

`summarise()` has grouped output by 'condition_rec'. You can override using the
`.groups` argument.

exp1_subj_filter_1 <- exp1_step1_goodsubj$Rec_Session_Id
exp1_item_filter_1 <- exp1_step1_item.err$target_rec

exp1_data_step1 <- exp1_rawdata.sub |>
  subset(condition_rec %in% exp1_info$freq_conditions) |>
  subset(Rec_Session_Id %in% exp1_subj_filter_1 & 
         !target_rec %in% exp1_item_filter_1)

Step 2: prime durations

In [6]:

exp1_summary.primeTime <- exp1_rawdata.sub %>% 
  summarise(meanPrimeTime = round(mean(primeTime), 2), 
            sdPrimeTime = round(sd(primeTime), 2))

exp1_primeTimeRangeSummary <- exp1_rawdata.sub %>% 
  group_by(primeTime) %>%
  mutate(range = ifelse(primeTime < exp1_info$prime_dur_lb, "below", 
                        ifelse(primeTime > exp1_info$prime_dur_ub, "above",
                               "in range"))) %>% 
  group_by(range) %>% tally() %>% ungroup() %>%
  mutate(range.percent = round((n*100)/nrow(exp1_rawdata.sub),2))

exp1_data_step2 <- exp1_data_step1  |>
  subset(primeTime >= exp1_info$prime_dur_lb & primeTime <= exp1_info$prime_dur_ub)

exp1_step2_subj_remain <- exp1_data_step2$Rec_Session_Id |>
  unique() |>
  length()

exp1_step2_trials_remain <- nrow(exp1_data_step2)

Step 3: RT analysis

In [7]:

# RT outliers 
exp1_data_step3 <- exp1_data_step2 |> 
  subset(RT >= exp1_info$rt_lb & RT <= exp1_info$rt_ub)

exp1_step3_subj_remain <- exp1_data_step3$Rec_Session_Id |>
  unique() |>
  length()

exp1_step3_trials_remain <- nrow(exp1_data_step3)

# error trial removal

exp1_data_step3b <- exp1_data_step3  |>
  subset(responseError == 0)

exp1_step3b_subj_remain <- exp1_data_step3b$Rec_Session_Id |>
  unique() |>
  length()

exp1_step3b_trials_remain <- nrow(exp1_data_step3b)

# remove subjects with less than 12 trials in at least one condition*primetype combination (half of the total number of items per combination)
rt_data_labels <- c("Rec_Session_Id", "condition_rec", "primetype_rec", "RT")

exp1_subj_filter_2 <- exp1_data_step3b[, rt_data_labels] |>
  aggregate(RT ~ ., FUN = length, drop = FALSE) |>
  subset(RT < 7, select = Rec_Session_Id) |>
  unique() |>
  unlist()

### we also want to sure that all subjects have all conditions; in case some subject had all the trials for a given condition lost down the road, they will be removed
exp1_subj_filter_conditions <- 
  exp1_data_step3b %>%
  group_by(Rec_Session_Id) %>% 
  distinct(condition_rec, primetype_rec) %>% 
  tally() %>% filter(n != 6) %>% pull(Rec_Session_Id)

exp1_data_final <- exp1_data_step3b |>
  subset(!(Rec_Session_Id %in% exp1_subj_filter_2) & !(Rec_Session_Id %in% exp1_subj_filter_conditions)) %>%
  mutate(condition_rec = as.factor(condition_rec), primetype_rec=as.factor(primetype_rec))

exp1_final_subj_remain <- exp1_data_final$Rec_Session_Id |>
  unique() |> 
  length()
  
exp1_final_trials_remain <- nrow(exp1_data_final)

Results

In [8]:

# error rates averages
### we also want to sure that all subjects have all conditions; in case some subject had all the trials for a given condition lost down the road, they will be removed. Crucially the trial calculations are made on the dataset *before* the trial error removal step
exp1_subj_filter_2_with.errors <- exp1_data_step3[, rt_data_labels] |>
  aggregate(RT ~ ., FUN = length, drop = FALSE) |>
  subset(RT < 7, select = Rec_Session_Id) |>
  unique() |>
  unlist()
# this step just makes sure that the same subjects will be removed from both datasets
exp1_subj_filter_2_with.errors <- union(exp1_subj_filter_2_with.errors, exp1_subj_filter_2)

### just making sure that all subjects have all conditions; in case some subject had all the trials for a given condition lost down the road, they will be removed
exp1_subj_filter_conditions_with.errors <- 
  exp1_data_step3 %>%
  group_by(Rec_Session_Id) %>% 
  distinct(condition_rec, primetype_rec) %>% 
  tally() %>% filter(n != 6) %>% pull(Rec_Session_Id)

exp1_data_final_with.errors <- exp1_data_step3 |>
  subset(!(Rec_Session_Id %in% exp1_subj_filter_2_with.errors) & 
           !(Rec_Session_Id %in% exp1_subj_filter_conditions_with.errors)) 

exp1_error.rates <- exp1_data_final_with.errors %>%
  mutate(primetype_rec = factor(primetype_rec, levels=c("unrelated", "related")),
         condition_rec = factor(condition_rec, levels=c("high", "low", "non-word"))) %>%
  group_by(condition_rec, primetype_rec, Rec_Session_Id) %>%
  summarise(error.percent=mean(responseError)*100)

`summarise()` has grouped output by 'condition_rec', 'primetype_rec'. You can
override using the `.groups` argument.

# RT averages
exp1_rt.avg_subj <- exp1_data_final %>% 
  group_by(Rec_Session_Id, condition_rec, primetype_rec) %>%
  summarise(meanRT = mean(RT))

`summarise()` has grouped output by 'Rec_Session_Id', 'condition_rec'. You can
override using the `.groups` argument.

exp1_rt_cor <- exp1_data_final %>% 
  group_by(Rec_Session_Id, condition_rec, primetype_rec) %>%
  dplyr::summarise(meanRT=mean(RT)) %>%
  select(Rec_Session_Id, condition_rec, primetype_rec, meanRT) %>% 
  pivot_wider(names_from='primetype_rec', values_from=c('meanRT')) %>%
  group_by(condition_rec) %>%
  summarise(cor=cor(unrelated, related))

`summarise()` has grouped output by 'Rec_Session_Id', 'condition_rec'. You can
override using the `.groups` argument.

# RT + error averages (by subject)
exp1_avgs_subj <- merge(exp1_rt.avg_subj, exp1_error.rates, by=c("Rec_Session_Id", "condition_rec", "primetype_rec"))

exp1_mop.err_across <- exp1_avgs_subj %>%
  group_by(condition_rec, primetype_rec) %>%
  summarise(gd.mean=mean(meanRT, na.rm=T), sd=sd(meanRT, na.rm=T), mean.error=mean(error.percent)) %>%
  pivot_wider(id_cols="condition_rec", names_from=primetype_rec, values_from=c(gd.mean, sd, mean.error)) %>%
  left_join(exp1_rt_cor, by="condition_rec")

`summarise()` has grouped output by 'condition_rec'. You can override using the
`.groups` argument.

#### ERROR PRIMING CALCULATIONS ####
# we will just run stats for this, no numerical calculations (e.g., priming)

### 2x2 ANOVA
exp1_errors_anova <- exp1_error.rates |>
  transform(Rec_Session_Id = factor(Rec_Session_Id)) |>
  subset(condition_rec %in% c("high", "low")) |>
  aov(error.percent ~ condition_rec * primetype_rec + Error(Rec_Session_Id/(condition_rec * primetype_rec)), data = _)

exp1_error_stats_main <- exp1_error.rates %>%
  mutate(primetype_rec = fct_relevel(primetype_rec, "unrelated", "related")) %>%
  group_by(condition_rec) %>%
  t_test(error.percent ~ primetype_rec, paired=T, detailed=T) %>% select(-3:-7) %>%
  rename(factor = "condition_rec", t='statistic')

### 2x2 GLMER
options(contrasts = c("contr.sum","contr.poly"))

exp1_error_2x2glmer <- glmer(responseError ~ condition_rec * primetype_rec + (1|Rec_Session_Id) + (1|target_rec), 
                             data = subset(exp1_data_final_with.errors, condition_rec %in% c("high", "low")), 
                             family = binomial, control=glmerControl(optimizer = "bobyqa", optCtrl=list(maxfun=1e6)))

Anova(exp1_error_2x2glmer)

Analysis of Deviance Table (Type II Wald chisquare tests)

Response: responseError
                               Chisq Df Pr(>Chisq)    
condition_rec                12.8737  1  0.0003332 ***
primetype_rec               421.3616  1  < 2.2e-16 ***
condition_rec:primetype_rec   4.9708  1  0.0257789 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

### 2x3 GLMER 
exp1_error_2x3glmer <- glmer(responseError ~ condition_rec * primetype_rec + (1|Rec_Session_Id) + (1|target_rec), data = exp1_data_final_with.errors, family = binomial, control=glmerControl(optimizer = "bobyqa", optCtrl=list(maxfun=1e6)))

Anova(exp1_error_2x3glmer, type = 3)

Analysis of Deviance Table (Type III Wald chisquare tests)

Response: responseError
                              Chisq Df Pr(>Chisq)    
(Intercept)                 3028.42  1  < 2.2e-16 ***
condition_rec                 14.61  2  0.0006722 ***
primetype_rec                253.46  1  < 2.2e-16 ***
condition_rec:primetype_rec  307.44  2  < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

exp1_error_emmeans_glmm <- emmeans(exp1_error_2x3glmer, pairwise ~ primetype_rec|condition_rec)

### MAIN EFFECTS
#### t-test
exp1_errors_stats_main <- exp1_error.rates %>% 
  group_by(condition_rec) %>%
  t_test(error.percent ~ primetype_rec, paired=T)

#### RT PRIMING CALCULATIONS ####

### MAIN EFFECTS

### 2x2 ANOVA
exp1_mop_anova <- exp1_rt.avg_subj |>
  transform(Rec_Session_Id = factor(Rec_Session_Id)) |>
  subset(condition_rec %in% c("high", "low")) |>
  aov(meanRT ~ condition_rec * primetype_rec + Error(Rec_Session_Id/(condition_rec * primetype_rec)), data = _)

### 2x2 GLMER 
exp1_mop_2x2glmer <- glmer(RT ~ condition_rec * primetype_rec + (1|Crowdsourcing_SubjId) + (1|target_rec), 
                        data = subset(exp1_data_final, condition_rec %in% c("high", "low")), 
                        family = Gamma(link="identity"), control=glmerControl(optimizer = "bobyqa", optCtrl=list(maxfun=1e6)))

Anova(exp1_mop_2x2glmer, type = 3)

Analysis of Deviance Table (Type III Wald chisquare tests)

Response: RT
                                 Chisq Df Pr(>Chisq)    
(Intercept)                 509615.850  1  < 2.2e-16 ***
condition_rec                  309.266  1  < 2.2e-16 ***
primetype_rec                 2275.770  1  < 2.2e-16 ***
condition_rec:primetype_rec     93.541  1  < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

### 2x3 GLMER 
exp1_mop_2x3glmer <- glmer(RT ~ condition_rec * primetype_rec + (1|Crowdsourcing_SubjId) + (1|target_rec), 
                        data = exp1_data_final %>% mutate(primetype_rec = factor(primetype_rec, levels=c("unrelated", "related"))),
                        family = Gamma(link="identity"), control=glmerControl(optimizer = "bobyqa", optCtrl=list(maxfun=1e6)))

Anova(exp1_mop_2x3glmer, type = 3)

Analysis of Deviance Table (Type III Wald chisquare tests)

Response: RT
                                Chisq Df Pr(>Chisq)    
(Intercept)                 2241361.8  1  < 2.2e-16 ***
condition_rec                  2521.4  2  < 2.2e-16 ***
primetype_rec                  1617.1  1  < 2.2e-16 ***
condition_rec:primetype_rec    1236.6  2  < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#estimates
exp1_emmeans_glmm <- emmeans(exp1_mop_2x3glmer, pairwise ~ primetype_rec|condition_rec) |> rbind() |> as_tibble() %>%
  select(-contrast, -df) %>% 
  pivot_wider(names_from=primetype_rec, values_from = c(emmean, SE, asymp.LCL, asymp.UCL)) %>%
  mutate(across(contains("emmean"), \(x) round(x)), 
         across(contains("asymp"), \(x) round(x)),
         across(contains("SE"), \(x) round(x, 2))) %>%
  relocate(c(emmean_unrelated, SE_unrelated, asymp.LCL_unrelated, asymp.UCL_unrelated), .before=emmean_related) %>%
    unite("CI_unrelated", asymp.LCL_unrelated:asymp.UCL_unrelated, sep=" ") %>%
  relocate(c(emmean_related, SE_related, asymp.LCL_related, asymp.UCL_related), .before=emmean_related) %>%
    unite("CI_related", asymp.LCL_related:asymp.UCL_related, sep=" ") %>%
    unite("CI_MOP", asymp.LCL_.:asymp.UCL_., sep=" ") %>%
  rename(MOP = "emmean_.", SE_MOP = "SE_.")

exp1_contrasts_glmm <- emmeans(exp1_mop_2x3glmer, pairwise ~ primetype_rec|condition_rec)$contrasts |> rbind() |> as_tibble() %>%
  select(condition_rec, z.ratio:p.value) %>% mutate(across(p.value, \(x) signif(x, 3)))

exp1_estimates_glmm <- exp1_emmeans_glmm %>%
  left_join(exp1_contrasts_glmm, by='condition_rec') %>%
  mutate(across(contains("CI"), \(x) paste0("[", x, "]"))) %>%
  rename(factor = 'condition_rec')

#### by subject
exp1_mop_subj <- exp1_rt.avg_subj %>% 
  pivot_wider(names_from=primetype_rec, values_from=meanRT) %>%
  mutate(priming = unrelated-related)

#### descriptive stats
exp1_gdavg_mop <- exp1_mop_subj %>%
  group_by(condition_rec) %>%
  summarise(MOP = mean(priming), se = sd(priming)/sqrt(n()), ci=(qt(0.975, n()-1)*se),
            sd=sd(priming), ES=round(MOP/sd, 2))

#### summary table
exp1_gdavg_mop_summary <- exp1_gdavg_mop %>%
  left_join(., exp1_mop.err_across, by='condition_rec') %>%
  mutate(across(c(2:5, 7:12), round), across(c(13), round, 2)) %>%
  mutate(ci.lb = paste0("[", MOP-ci), ci.ub = paste0(MOP+ci, "]")) %>%
  unite("CI", ci.lb:ci.ub, sep = " ") %>% select(-ci, -se) %>% 
  rename(factor = "condition_rec") %>%
  relocate(gd.mean_related:cor, .before=MOP) %>% 
  relocate(gd.mean_related, .after=gd.mean_unrelated) %>%
  relocate(CI, .after=MOP)

Warning: There was 1 warning in `mutate()`.
ℹ In argument: `across(c(13), round, 2)`.
Caused by warning:
! The `...` argument of `across()` is deprecated as of dplyr 1.1.0.
Supply arguments directly to `.fns` through an anonymous function instead.

  # Previously
  across(a:b, mean, na.rm = TRUE)

  # Now
  across(a:b, \(x) mean(x, na.rm = TRUE))

#### t-test
exp1_rt_stats_main <- exp1_rt.avg_subj %>%
  mutate(primetype_rec = fct_relevel(primetype_rec, "unrelated", "related")) %>%
  group_by(condition_rec) %>%
  t_test(meanRT ~ primetype_rec, paired=T) %>% select(-2:-6) %>%
  rename(factor = "condition_rec", t='statistic') %>%
  mutate_if(is.numeric, list(~as.character(signif(., 3))))

### INTERACTION EFFECT

#### descriptive stats
exp1_gdavg_fae <- exp1_mop_subj %>%
  select(-unrelated, -related) %>%
  pivot_wider(names_from="condition_rec", values_from="priming") %>%
  mutate(interaction = low-high) %>% ungroup() %>%
  summarise(factor = "frequency:primetype", 
            sd=sd(interaction), se = sd/sqrt(n()), ci = (qt(0.975, n()-1)*se)) %>%
  mutate(across(c(2:4), round))

##### summary table
exp1_gdavg_fae_summary <- exp1_mop_subj %>% 
  select(-unrelated, -related) %>% ungroup() %>%
  pivot_wider(names_from="condition_rec", values_from="priming") %>%
  summarise(mean_high = round(mean(high)), mean_low = round(mean(low)), cor=round(cor(high, low), 3)) %>%
  bind_cols(., exp1_gdavg_fae) %>%
  mutate(MOP = mean_low-mean_high, ES = round(MOP/sd, 2), ci.lb = paste0("[", MOP-ci), ci.ub = paste0(MOP+ci, "]")) %>%
  unite("CI", ci.lb:ci.ub, sep = " ") %>% select(-ci, -se) %>% 
  relocate(CI, .after=MOP) %>%
  relocate(factor, .before=mean_high)
  
#### t-test (from the ANOVA)
exp1_rt_stats_interaction <- exp1_mop_subj %>% ungroup() %>%
  filter(condition_rec != 'non-word') %>%
  mutate(condition_rec = fct_relevel(condition_rec, "low", "high")) %>%
  t_test(priming ~ condition_rec, paired=TRUE) %>% 
  mutate_if(is.numeric, list(~as.character(signif(., 3)))) %>%
  select(-1, -4:-5) %>%
  unite("factor", group1:group2) %>%
  rename(t='statistic') %>%
  mutate(factor = "frequency:primetype")

## estimates from the GLMM model
exp1_emmeans_glmm_fae <- confint(contrast(emmeans(exp1_mop_2x3glmer, ~ primetype_rec*condition_rec), interaction = c("pairwise", "pairwise"), adjust='bonferroni')) %>% slice(1) %>% select(estimate:SE, asymp.LCL, asymp.UCL) %>% 
  mutate(across(c(1, 3:4), \(x) x*-1),
         across(c(1:4), \(x) round(x,2))) %>%
  mutate(factor = "frequency*primetype") %>% rename(MOP = "estimate", SE_MOP = "SE") %>% unite("CI_MOP", asymp.UCL:asymp.LCL, sep=" ") %>%
  relocate(CI_MOP, .after=SE_MOP) %>% relocate(factor, .before=MOP) %>% mutate(CI_MOP = paste0("[", CI_MOP, "]"))

exp1_z_glmm_fae <- contrast(emmeans(exp1_mop_2x3glmer, ~ primetype_rec*condition_rec), interaction = c("pairwise", "pairwise"), adjust='bonferroni') |>
  rbind() |> as_tibble() %>% slice(1) %>% select(z.ratio:p.value) %>% 
  mutate(across(c(1), \(x) x*-1)) %>%
  mutate(factor = "frequency*primetype") %>% relocate(factor, .before=z.ratio) 

exp1_estimates_glmm_fae <- left_join(exp1_emmeans_glmm_fae, exp1_z_glmm_fae, by='factor')

Stats summary

In [9]:

In [10]:

exp1_summary.results_mop <- merge(exp1_gdavg_mop_summary, exp1_rt_stats_main, by='factor')
exp1_summary.results_fae <- merge(exp1_gdavg_fae_summary, exp1_rt_stats_interaction, by='factor') |>
  select(-mean_high, -mean_low)

exp1_summary.results <- bind_rows(exp1_summary.results_mop, exp1_summary.results_fae)
  
exp1_summary.results %>%
  relocate(c("sd_unrelated", "mean.error_unrelated"), .before=gd.mean_related) %>%
  gt() %>%
  cols_label(
    CI = "95% CI",
    contains("mean") ~ "mean",
    contains("sd") ~ "SD", 
    contains("error") ~ "Error (%)"
  ) %>%
  tab_spanner(
    label = "unrelated RT",
    columns = c(2:4)
  ) %>%
  tab_spanner(
    label = "repetition RT",
    columns = c(5:7)
  ) %>%
  tab_spanner(
    label = 'priming effects',
    columns = c(9:12)
  ) %>%
  tab_spanner(
    label = md("_t_-test"),
    columns = c(13:15)
  ) %>%
  cols_label(
    sd = md("SD~p~")
  ) %>%
  cols_label(
    t = md("_t_"),
    p = md("_p_"),
  ) %>%
   sub_missing(
    missing_text = " "
  )

Experiment 1. Summary of the word priming results within the ANOVA model. Legend. MOP: magnitude of priming.

factor	unrelated RT			repetition RT			cor	priming effects				t-test
	mean	SD	Error (%)	mean	SD	Error (%)		MOP	95% CI	SDp	ES	t	df	p
high	573	83	3	555	85	2	0.860	18	[16 20]	45	0.41	19.7	2340	2.88e-80
low	605	88	6	577	88	3	0.850	28	[26 30]	49	0.58	27.8	2340	1.52e-147
non-word	623	103	4	625	103	4	0.910	-2	[-4 0]	43	-0.05	-2.33	2340	0.0197
frequency:primetype							0.029	10	[7 13]	66	0.15	7.24	2340	5.86e-13

In [11]:

In [12]:

exp1_summary.results_glmm <- bind_rows(exp1_estimates_glmm, exp1_estimates_glmm_fae)
  
exp1_summary.results_glmm %>%
  mutate(z.ratio = round(z.ratio, 2),
         p.value = "<.0001") %>%
  gt() %>%
  cols_label(
    contains("CI") ~ "95% CI",
    contains("emmean") ~ "mean",
    contains("SE") ~ "SE"
  ) %>%
  tab_spanner(
    label = "unrelated RT",
    columns = c(2:4)
  ) %>%
  tab_spanner(
    label = "repetition RT",
    columns = c(5:7)
  ) %>%
  tab_spanner(
    label = 'priming effects',
    columns = c(8:10)
  ) %>%
   sub_missing(
    missing_text = " "
  )

Experiment 1. Summary of the word priming results based on the 2x3 GLMM model, with the same fixed and random effects as the 2x2 GLMM model, but including the non-word condition. Legend. MOP: magnitude of priming.

factor	unrelated RT			repetition RT			priming effects			z.ratio	p.value
	mean	SE	95% CI	mean	SE	95% CI	MOP	SE	95% CI
high	604	0.78	[602 606]	585	0.77	[583 587]	19.00	0.59	[17 21]	32.21	<.0001
low	634	0.62	[633 636]	606	0.65	[604 608]	28.00	0.65	[27 30]	43.51	<.0001
non-word	645	0.82	[643 647]	647	0.79	[645 649]	-2.00	0.57	[-4 0]	-3.42	<.0001
frequency*primetype							9.42	0.77	[7.59 11.26]	12.27	<.0001