---
title: Familiar sequences are processed faster than unfamiliar sequences, even when
  they do not match the count-list
output:
  pdf_document: default
  html_document:
    df_print: paged
author: Devlin et al. (2024)
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, fig.align = 'center', out.width = '70%')
```

```{r load packages, include=FALSE, echo=FALSE}
library(tidyverse)
library(readxl)
library(sirt)
library(moments)
library(broom)
library(lme4)
library(ggcorrplot)
library(grid)
library(gridExtra)

rm(list = ls())

# Open the RProj file "setWD" to set the working directory to this folder before running the script.
```

```{r load data, include=FALSE, echo=FALSE}
order_df <- read_excel("data.xlsx", sheet = "order")
CJ_df <- read_excel("data.xlsx", sheet = "CJ")

```

# Analysis

## Data quality and accuracy

60 participants completed a comparative judgement task, an order verification task, and an arithmetic fluency task.

We first checked accuracy on the order verification task to establish whether participants understood the task instructions and completed the task correctly.

An example of the order verification task data is shown below:

```{r show dataset, echo=FALSE}
knitr::kable(head(order_df, 5), "simple", digits = 2, caption = "First five rows of order dataset")
```

We first checked participants overall accuracy on the order verification task:

```{r calculate accuracy, echo=TRUE}
percent_overall <- order_df %>%
  group_by(ppt) %>% 
  count(corr_resp) %>%
  mutate(percent = n / sum(n)) %>%
  filter(corr_resp == 1) %>%
  select(ppt, percent_correct = percent)
```

```{r check accuracy, echo=TRUE}
macc <- mean(percent_overall$percent_correct)
sdacc <- sd(percent_overall$percent_correct)

lowerbound <- macc - (sdacc * 3)

below3 <- percent_overall %>%
  filter(percent_correct < lowerbound)

```

This revealed that two participants displayed unusually low accuracy (i.e., 3 SDs below mean):

```{r acc 3SDs below mean, echo=FALSE}
knitr::kable(below3, "simple", digits = 2, caption = "Accuracy 3 SDs below the mean")
```

Furthermore, we also observed two other participants to display below chance accuracy on certain sequence types: 

```{r accuracy sequence types, echo=TRUE}
percent_type <- order_df %>%
  group_by(ppt, type) %>% 
  count(corr_resp) %>%
  mutate(percent = n / sum(n)) %>%
  filter(corr_resp == 1) %>%
  select(ppt, percent_correct = percent)

below_chance <- percent_type %>%
  filter(percent_correct < 0.5)

```

```{r acc below chance, echo=FALSE}
knitr::kable(below_chance, "simple", digits = 2, caption = "Accuracy below chance")
```

Therefore, data from these participants were removed from all tasks - resulting in 56 participants being included in the final analyses:

```{r remove ppts}
order_df <- order_df %>% 
  filter(!ppt %in% c("220", "215", "610", "617"))

CJ_df <- CJ_df %>% 
  filter(!judgeName %in% c("220", "215", "610", "617"))
```

## Comparative judgement analysis

```{r CJ analysis prepare data, echo=FALSE}
decisions <- CJ_df
```

To produce a ranking and relative scoring of sequences, the raw judgement data (i.e., comparisons) from all participants were fitted to the Bradley-Terry (1952) model. 

```{r run CJ anlysis, message=FALSE, echo=FALSE, results='hide', warning=FALSE}
# N.B. CJ code provided by Dr Ian Jones (Department of Mathematics Education, Loughborough University, UK).
# prepare data for fitting to the BT model
players <- unique(c(decisions$chosen,decisions$notChosen))
df <- data.frame(id1=decisions$chosen, id2=decisions$notChosen,  result=1, judge=decisions$judgeName)
df$id1 <- factor(df$id1,levels=players)
df$id2 <- factor(df$id2,levels=players)

# fit to BT model and get scale values and SSR
mod1 <- sirt::btm(df[,1:3]   , judge = df$judge, maxit=400 , fix.eta=0 , ignore.ties=TRUE )
summary(mod1) #scale values are the theta column in 
SSR <- mod1$mle.rel #Scale Separation Reliability: a bit like Cronbach's alpha

# estimate inter-rater reliability using the split-halves technique
iterations = 1000 #100 is usually fine, but change if you wish
ir = data.frame(iteration = rep(0,iterations), pearson = rep(0,iterations), SSR1 = rep(0,iterations), SSR2 = rep(0,iterations), num_cand = rep(0,iterations)) #stores the result of each iteration
ip = as.data.frame(matrix(data=0,nrow=nrow(mod1$effects),ncol=2*iterations)) # store all theta values in case of interest
colnames(ip) = as.character(c(1:(2*iterations)))
noj = length(unique(decisions$judgeName)) # number of judges
judges = data.frame(row = 1:noj, judgeName = unique(decisions$judgeName), rand = rep(0,noj), group = rep(0,noj)) #list of unique judge names

# loop 'iterations' times
for(i in 1:iterations) 
{ 
  ir[i,1] = i # fill iteration column
  
  # randomise judges to two groups
  judges$rand = runif(noj, -1, 1)  
  judges <- judges[order(judges$rand),]
  judges$row = 1:noj
  judges$group = 2
  judges$group[1:(noj/2)] <-1
  groups = merge(decisions,judges) # now split judgements corresponding to judge groups
  g1 = groups[groups$group==1,]
  g2 = groups[groups$group==2,]
  
  # calculate scale values for each group
  players1 <- unique(c(g1$chosen,g1$notChosen))
  players2 <- unique(c(g2$chosen,g2$notChosen))
  df1 <- data.frame(id1=g1$chosen, id2=g1$notChosen, result=1)
  df1$id1 <- factor(df1$id1,levels=players)
  df1$id2 <- factor(df1$id2,levels=players)
  modg1 <- sirt::btm(df1 , maxit=400 , fix.eta=0 , ignore.ties=TRUE )
  
  df2 <- data.frame(id1=g2$chosen, id2=g2$notChosen, result=1)
  df2$id1 <- factor(df2$id1,levels=players)
  df2$id2 <- factor(df2$id2,levels=players)
  modg2 <- sirt::btm(df2 , maxit=400 , fix.eta=0 , ignore.ties=TRUE )
  
  modg12 = merge(modg1$effects, modg2$effects, by="individual") # bodge in case different number of candidates in each rank
  
  
  # calculate Pearson correlation and store results
  ir[i,2] =  cor(modg12$theta.x,modg12$theta.y, method="pearson") # fill correlation column
  ir[i,3] =  modg1$mle.rel
  ir[i,4] =  modg2$mle.rel 
  ir[i,5] = nrow(modg12)
  ip[,2*i-1] = c(modg12$theta.x, rep(0,nrow(mod1$effects)-nrow(modg12)))
  ip[,2*i] = c(modg12$theta.y, rep(0,nrow(mod1$effects)-nrow(modg12)))
  
  
  # display iteration number so user knows it's not crashed
  print(i)
} # end of loop

# take median value
interrater <- median(ir$pearson)

# select columns from output
fam_scores <- mod1$effects %>%
  select(seq = individual, fam_score = theta)
```

For ease of interpretation, we then re-scaled these scores between 0 and 100:

```{r scale familiarity scores between 0 and 100}
# Scale theta scores to between 0 and 100
scale <- function(x) (x-min(x))/(max(x) - min(x)) * 100
fam_scores$scaled_score <- scale(fam_scores$fam_score)
```

```{r calculate SSR and IRR, echo = FALSE, message = FALSE, result = FALSE}
SSR <- round(SSR, 2)
IRR <- round(interrater, 2)
```

To establish the reliability of these judgements, we considered two reliability measures common in the comparative judgement literature. First, we calculated scale separation reliability – a measure of internal consistency analogous to Cronbach’s alpha. This returned a score of .997, indicating very good internal consistency. 

Second, we considered inter-rater reliability – a reliability measure calculated by repeatedly splitting the judges (i.e., participants) into two groups and then calculating the correlation between their produced scores. 

Following this procedure with 1000 repetitions, we observed a median Pearson product-moment correlation coefficient of .99, suggesting very good inter-rater reliability. Overall, these measures suggest a high degree of agreement across participants in their judgements. 

The familiarity scores generated by the comparative judgement task can be seen below:

```{r add distance and direction to famscores, echo=FALSE}
type <- read_excel("data.xlsx", sheet = "seqDetails")

fam_scores <- fam_scores %>% 
  left_join(type, by = 'seq')
```

```{r familiarity scores plot, message=FALSE, echo=FALSE, out.width="95%", warning=FALSE}
ggplot(fam_scores, aes(scaled_score, reorder(seq, scaled_score), color = distance)) +
  geom_point(size = 4, shape = 16, aes(fill = direction), show.legend = FALSE) +
  geom_segment(aes(xend = 0, yend = seq, linetype = direction), size = 1.5) +
  geom_text(aes(label = seq), color = "black", size = 4, hjust = -.3, vjust = .34, fontface = "bold") +
  xlim(0,105) +
  labs(title = element_blank(), x = "Familiarity scores (0-100)", y = " ") +
  theme_classic() +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        legend.position = c(0.9, 0.3)) +
  scale_color_manual(values = c("deepskyblue1", "#FF6A6A"), name = "", labels = c("Consecutive", "Non-consecutive")) +
  scale_linetype_manual(values = c(1,6), name = " ", labels = c("Ascending", "Descending"))
```

```{r means, echo=FALSE, results=FALSE}
fam_scores %>%
  select(seq, scaled_score, size)
```

```{r means for different sequence types}
fam_scores %>% 
  group_by(distance) %>% 
  summarise(mean = mean(scaled_score))

fam_scores %>% 
  group_by(direction) %>% 
  summarise(mean = mean(scaled_score))
```

Consistent with our expectations, consecutive sequences (mean familiarity score = 55.65) were judged as more familiar than non-consecutive sequences (mean = 23.08). Similarly, ascending sequences (mean = 47.20) were judged as more familiar than descending sequences (mean = 29.36). 

Furthermore, we note that 1-2-3 was judged as the most familiar sequence (familiarity score = 100), which is consistent with suggestions this sequence is processed atypically fast due to its high familiarity (i.e., Sella et al., 2020). 

We also observed 2-4-6 to be the most familiar non-consecutive sequence (familiarity score = 51.44). In fact, 2-4-6 was judged as being more familiar than over half (8 out of 14) of the consecutive sequences. 

Similarly, 3-6-9 (familiarity score = 43.25) and 4-6-8 (familiarity score = 40.72) had relatively high familiarity scores despite being non-consecutive. This is consistent with suggestions that sequences that are a part of highly practiced multiplication tables are likely to be highly familiar (i.e., Vos et al., 2021).

## Order verification analysis

Each participant completed 600 order verification trials, meaning 33,600 responses were recorded across 56 participants. 

```{r calculate mean accuracy}
percent_overall <- order_df %>%
  group_by(ppt) %>% 
  count(corr_resp) %>%
  mutate(percent = n / sum(n)) %>%
  filter(corr_resp == 1) %>%
  select(ppt, percent_correct = percent)

mean_accuracy <- round(mean(percent_overall$percent_correct * 100),2)
```

Mean accuracy across all participants was 93.55%.

Because we were only interested in correct responses for ordered sequences, we first removed all non-ordered trials (n = 16,800; 50%) and all incorrect responses (n = 1,027; 6.11%). 

```{r remove incorrect and non-ordered responses}
# remove non-ordered sequences
order_df <- order_df %>% 
  filter(order == 1)

# remove incorrect responses
order_df <- order_df %>% 
  filter(corr_resp == 1)
```

Consistent with previous studies (e.g., Dubinkina et al., 2021) we also removed responses that were shorter than 200 ms (n = 8; 0.05%) or longer than the overall mean response time plus three standard deviations (mean = 1432 ms, SD = 1598 ms) (n = 225; 1.43%) as these were considered unlikely to be genuine responses:

```{r remove below 200 ms or above mean + 3 SDs}
# remove RTs below 200 ms (consistent with Dubinkina21)
order_df <- order_df %>% 
  filter(RT > 0.2)

# remove RTs greater than the mean + 3 SDs (consistent with Dubinkina21)
order_df <- order_df %>% 
  filter(RT < mean(RT) + 3*(sd(RT)))
```

Following this trimming procedure, the distribution appeared somewhat positively skewed (skewness = 2.4) and leptokurtic (kurtosis = 13.17):

```{r plot distribution}
ggplot(order_df, aes(x=RT)) +
  geom_histogram(fill = '#6495ED', alpha = 1, bins = 40) +
  labs(title = "Distribution of RTs", x = "RTs") +
  theme_classic()
```

```{r calculate skewness and kurtosis, echo=FALSE, results='hide'}
skewval <- skewness(order_df$RT)
kurtval <- (kurtosis(order_df$RT)) + 3

round(skewval, 2)
round(kurtval, 2)
```

To account for this, unless otherwise specified, we considered the median response times for each participant. Means and standard deviations were then calculated for each sequence type:

```{r means and SDs, message = FALSE}
# descriptives: means and sds for AC, DC, ANC, DNC
sum_df <- order_df %>% 
  group_by(ppt, consec, direction) %>% 
  summarise(median = median(RT) * 1000) %>% 
  ungroup()

descriptives <- sum_df %>% 
  group_by(consec, direction) %>% 
  summarise(mean = round(mean(median),0), sd = round(sd(median),0)) %>% 
  ungroup()

```

```{r means and SDs table, echo=FALSE}
knitr::kable(descriptives, "simple", digits = 2, caption = "Means and standard deviations for all direction/distance combinations (ms)")
```

### Relationship between response times and familiarity scores

Next, we considered the relationship between response times for each sequence on the order verification task and their associated familiarity scores generated by the comparative judgement task. 

```{r prepare dataframes for analyses, echo=FALSE, message=FALSE}
# create a high/low familiarity variable (split-half)
fam_splithalf <- fam_scores %>% 
  arrange(desc(scaled_score)) %>% 
  mutate(rank = row_number()) %>% 
  mutate(familiarity = case_when(rank <= 15 ~ "high",
                                 rank > 15 ~ "low")) %>% 
  select(seq, familiarity)

# calculate median RTs for each sequence for each participant
medRTs <- order_df %>% 
  group_by(ppt, seq) %>% 
  summarise(medRT = median(RT) * 1000) %>% 
  ungroup()

# take the mean median for each sequence
medRTs <- medRTs %>% 
  group_by(seq) %>% 
  summarise(meanRT = mean(medRT))

# join familiarity scores to RT data
medRTs <- medRTs %>%
  left_join(fam_scores, by = "seq") 

#medRTs <- medRTs %>%
#  left_join(type, by = "seq")

# add familiarity variable
medRTs <- medRTs %>% 
  left_join(fam_splithalf, by = "seq")
```

```{r corr, message = FALSE}
RTcorrFam <- cor.test(medRTs$meanRT, medRTs$scaled_score)
tidy(RTcorrFam)
```

We first conducted a Pearson’s product-moment correlation which revealed a significant negative correlation, r(28) = -0.82, p < .001, suggesting higher familiarity scores were associated with faster response times:

```{r famXRT plot, message=FALSE, echo=FALSE, out.width="95%", warning=FALSE}
linear_model <- lm(meanRT ~ scaled_score, data = medRTs)
r_squared <- summary(linear_model)$r.squared

x_position <- 6
y_position <- 700


ggplot(medRTs, aes(scaled_score, meanRT)) +
  geom_point(aes(color = distance, shape = direction), size = 2) +
  geom_smooth(method = 'lm', se = FALSE, color = "black") +
  theme_classic() +
  theme(legend.position = c(0.91, 0.8)) +
  scale_color_manual(values = c("deepskyblue1", "#FF6A6A"), name = element_blank(), labels = c("Consecutive", "Non-consecutive")) +
  scale_shape_manual(values = c(16,15), name = element_blank(), labels = c("Acending", "Descending")) +
  labs(title = element_blank(), x = "Familiarity score (0-100)", y = "Response time (ms)") +
  geom_text(data = subset(medRTs, seq == "1 2 3" | seq == "2 4 6" | seq == "3 6 9" | seq == "3 2 1"),
            aes(label = seq), nudge_x = 5, nudge_y = 25, size = 4) +
  geom_text(aes(x = x_position, y = y_position, label = paste(expression(paste("R"^2,"= 0.67")))),
            col = "blue", size = 6, fontface = "bold", inherit.aes = FALSE, parse = TRUE)
```


However, although we observed a significant correlation between response times and familiarity scores, it could be argued this may simply be an artefact of ascending and consecutive sequences generally being more familiar than descending and non-consecutive sequences. 

Consequently, rather than demonstrating a familiarity effect, this association may instead simply reflect the known findings that consecutive sequences are generally processed faster than non-consecutive sequences and ascending sequences are generally processed faster than descending sequences. 

Therefore, we applied a linear mixed-effects modelling approach to establish whether familiarity scores explained any unique variance in response times beyond the effects of distance and direction alone.

Linear mixed-effects models have several advantages over more conventional approaches such as using ANOVAs or multiple regressions. For instance, not only are they compatible with both continuous and categorical explanatory variables, but they also allow taking into account variability across both participants and items simultaneously (Brown, 2021). 

Therefore, using the lme4 package in R (Bates et al., 2015), we examined the fixed effects of distance (consecutive/non-consecutive), direction (ascending/descending), and familiarity (0–100) while accounting for the random effects of participants and sequences.

N.B., because mixed effects models use unaggregrated data, median RTs are not used here.

```{r prepare data for mixed effects modelling, message=FALSE, echo=FALSE}
# unaggregated data
df_unaggregated <- order_df %>% 
  select(ppt, seq, RT, distance = consec, direction) %>% 
  mutate(RT = round(RT*1000, 0))

# recode (0 = C, 1 = NC; 0 = Asc, 1 = Dsc)
df_unaggregated <- df_unaggregated %>%
  mutate(distance = factor(ifelse(distance == "C", 0, 1)),
         direction = factor(ifelse(direction == "asc", 0, 1)))

fam_scores <- fam_scores %>% 
  select(-distance, -direction)

# add familiarity scores
df_unaggregated <- df_unaggregated %>% 
  left_join(fam_scores, by = c("seq"))




```


```{r run mixed effects models, message=FALSE}
# create lme4 models
mdl_null <- lmer(RT ~ 1 + (1 | ppt) + (1 | seq), 
                 data = df_unaggregated)

mdl_dist <- lmer(RT ~ distance + (1 | ppt) + (1 | seq), 
                 data = df_unaggregated)

mdl_dir <- lmer(RT ~ direction + (1 | ppt) + (1 | seq), 
                data = df_unaggregated)

mdl_fam <- lmer(RT ~ scaled_score + (1 | ppt) + (1 | seq), 
                data = df_unaggregated)

mdl_dist_dir <- lmer(RT ~ distance + direction + (1 | ppt) + (1 | seq), 
                     data = df_unaggregated)

mdl_dist_dir_fam <- lmer(RT ~ distance + direction + scaled_score + (1 | ppt) + (1 | seq), 
                         data = df_unaggregated)

# model comparisons
anova1 <- anova(mdl_null, mdl_dist_dir)
anova2 <- anova(mdl_dist_dir, mdl_dist_dir_fam)
anova3 <- anova(mdl_dist_dir_fam, mdl_fam)
```

We assessed model fit using likelihood-ratio tests. This revealed a model including distance and direction to better fit the data than a null model, $\chi^2$(2) = 8.88, p = .012:

```{r comparison1, echo=FALSE}
knitr::kable(anova1, "simple", digits = 2, caption = "Model comparison")
```

Furthermore, adding familiarity as a predictor further improved model fit, $\chi^2$(1) = 40, p < .001; thus, suggesting familiarity scores explained unique variance in response times beyond the effects of distance and direction alone:

```{r comparison2, echo=FALSE}
knitr::kable(anova2, "simple", digits = 2, caption = "Model comparison")
```

Moreover, when comparing the model incorporating distance, direction, and familiarity with a model including only familiarity, we found that the former still explained more variance ($\chi^2$(2) = 17.06, p < .001).

However, the BIC for the model with only familiarity (BIC=249442) was smaller than the model with distance, direction, and familiarity (BIC=249444), suggesting that the model with only familiarity was the most parsimonious explanation of response times in the number order verification task:

```{r comparison3, echo=FALSE}
knitr::kable(anova3, "simple", digits = 2, caption = "Model comparison")
```


As an exploratory analysis, we also re-ran these models incorporating numerical size 
as a predictor. Here, we calculated the numerical size as the sum of the items within each 
sequence (e.g., 1-2-3 has a numerical size of 6, or an average numerical size of 2):

```{r LME models including numerical size, message=FALSE}
# create lme4 models including size as a predictor
mdl_size <- lmer(RT ~ size + (1 | ppt) + (1 | seq), 
                data = df_unaggregated)

mdl_dist_dir_size <- lmer(RT ~ distance + direction + size + (1 | ppt) + (1 | seq), 
                     data = df_unaggregated)

mdl_dist_dir_size_fam <- lmer(RT ~ distance + direction + size + scaled_score + (1 | ppt) + (1 | seq), 
                         data = df_unaggregated)


# model comparisons
anova1 <- anova(mdl_null, mdl_dist_dir_size)
anova2 <- anova(mdl_null, mdl_dist_dir_size_fam)
anova3 <- anova(mdl_dist_dir_size_fam, mdl_fam)
```

We again assessed model fit using likelihood-ratio tests. This revealed a model including distance, direction, and numerical size to better fit the data than a null model, $\chi^2$(3) = 34.47, p < .001:

```{r comparison1 new, echo=FALSE}
knitr::kable(anova1, "simple", digits = 2, caption = "Model comparison")
```

Furthermore, adding familiarity as a predictor further improved model fit, $\chi^2$(4) = 55.05, p < .001; thus, suggesting familiarity scores explained unique variance in response times beyond the effects of distance, direction, and numerical size alone:

```{r comparison2 new, echo=FALSE}
knitr::kable(anova2, "simple", digits = 2, caption = "Model comparison")
```

Moreover, when comparing the model incorporating distance, direction, numerical size, and familiarity with a model including only familiarity, we found that the former still explained more variance ($\chi^2$(3) = 23.23, p < .001).

However, as in the previous analysis, the BIC for the familiarity-only model (BIC=249442) was smaller than the model with distance, direction, numerical size, and familiarity (BIC=249448), suggesting that the model with only familiarity was the most parsimonious explanation of response times in the number order verification task:

```{r comparison3 new, echo=FALSE}
knitr::kable(anova3, "simple", digits = 2, caption = "Model comparison")
```


### Reverse distance effect vs familiarity effect

```{r prepare dataframes, echo=FALSE, message=FALSE}
new_df <- order_df %>% 
  left_join(fam_splithalf, by = "seq")

# create RDE df
RDE_df <- new_df %>% 
  group_by(ppt, consec) %>% 
  summarise(median = median(RT) * 1000) %>% 
  ungroup()

RDE_wide <- RDE_df %>% 
  pivot_wider(names_from = consec, values_from = median)

RDE_df <- RDE_wide %>% 
  group_by(ppt) %>% 
  mutate(RDE = NC-C)

# creating FE variable
FE_df <- new_df %>% 
  group_by(ppt, familiarity) %>%
  summarise(medRT = median(RT)*1000)

FE_wide <- FE_df %>% 
  pivot_wider(names_from = familiarity, values_from = medRT)

FE_df <- FE_wide %>% 
  group_by(ppt) %>% 
  mutate(FE = `low` - `high`)

# joining RDE and FE dfs
RDE_FE <- RDE_df %>% 
  left_join(FE_df, by = "ppt") %>% 
  select(ppt, RDE, FE)

RDE_FE <- RDE_FE %>% 
  pivot_longer(cols = c("RDE", "FE"),
               names_to = "effect",
               values_to = "size")

RDE_FE$effect <- ifelse(RDE_FE$effect == "RDE", "Reverse distance effect", "Familiarity effect")
```

Next, we considered the reverse distance effect and its stability. We first calculated skewness and kurtosis individually for consecutive, non-consecutive, high familiarity, and low familiarity sequences. Although skewness was within acceptable bounds (skewness = 0.64-0.79), all sequence types were leptokurtic (kurtosis = 5.38-5.66). Therefore, we also report non-parametric results for all analyses in this section.

We then checked whether the reverse distance effect across both ascending and descending sequences was present at the group level.

```{r RDE plot, message = FALSE, out.width = "95%", echo=FALSE}
ggplot(RDE_df, aes(RDE, reorder(as.factor(ppt), RDE))) +
  geom_point(size=2, shape=21, fill="red", alpha=0.6) +
  geom_vline (xintercept=0, linetype='dashed', color='grey1', alpha=0.3) +
  coord_cartesian(xlim =c(-2000, 2000)) +
  labs(title="RDExPpt", x="Reverse distance effect (ms)", y=element_blank()) +
  theme_classic()
```

A paired-samples t-test revealed responses were indeed faster for consecutive sequences (mean = 1056 ms, SD = 312) than for non-consecutive sequences (mean = 1144 ms, SD = 373), t(55) = 3.99, p < .001 [V = 331, p < .001]:

```{r RDE stability, message = FALSE}
# t-test for RDE
RDE_ttest <- t.test(RDE_df$C, RDE_df$NC, paired = TRUE)
tidy(RDE_ttest)

```

```{r RDE stability non-para, message = FALSE}
# Wilcoxon Signed rank test for RDE
RDE_WCtest <- wilcox.test(RDE_df$C, RDE_df$NC, paired = TRUE)
tidy(RDE_WCtest)
```

```{r RDE stability2, echo = FALSE, message = FALSE, results='hide'}
round(mean(RDE_df$C),0)
round(sd(RDE_df$C),0)
round(mean(RDE_df$NC),0)
round(sd(RDE_df$NC),0)
```

Following this, we considered the stability of the reverse distance effect at the individual level. Therefore, we calculated the raw differences in response times between consecutive and non-consecutive sequences for each participant (mean = 88 ms, SD = 165):

```{r RDE and FE plot, message = FALSE, out.width = "95%", echo=FALSE}
ggplot(RDE_FE, aes(size, reorder(as.factor(ppt), size), fill = effect)) +
  geom_point(size=2, shape=21, alpha=0.6) + 
  geom_vline (xintercept=0, linetype='dashed', color='grey1', alpha=0.4) +
  coord_cartesian(xlim =c(-1200, 1200)) +
  labs(title= element_blank(), x="Difference in ms", y="Participants") +
  scale_fill_manual("Effect", values = c("green", "red")) +
  scale_x_continuous(breaks = c(-1000,-500, 0, 500, 1000)) +
  theme_classic() +
  theme(legend.position = c(.85, .2)) +
  theme(axis.text.y = element_blank())
```

We then considered the possibility of a familiarity effect. To calculate this, we used the familiarity scores from the comparative judgement task; the top half of sequences were coded as familiar, and the bottom half were coded as unfamiliar. 

```{r FE analyses, message = FALSE}
# t-test for familiarity effect
FE_ttest <- t.test(FE_df$low, FE_df$high, paired = TRUE)
tidy(FE_ttest)
```

```{r FE analyses non-para, message = FALSE}
# Wilcoxon Signed rank test for FE
FE_WCtest <- wilcox.test(FE_df$low, FE_df$high, paired = TRUE)
tidy(FE_WCtest)
```

```{r FE analyses2, message = FALSE, echo=FALSE, results='hide'}
round(mean(FE_df$high),0)
round(sd(FE_df$high),0)
round(mean(FE_df$low),0)
round(sd(FE_df$low),0)
```

A paired-samples t-test revealed response times were faster for familiar sequences (mean = 985 ms, SD = 285) than for unfamiliar sequences (mean = 1243 ms, SD = 428), t(55) = 9.07, p < .001 [V = 1596, p < .001]; thus, indicating a familiarity effect was present at the group level. 

We then considered the stability of this familiarity effect at the individual level. In contrast to the reverse distance effect, we observed that every participant showed a familiarity effect larger than 0 ms (mean = 258 ms, SD = 213). 

```{r RDE vs FE, message = FALSE}
# difference between RDE and FE
RDE_FE1 <- RDE_df %>% 
  left_join(FE_df, by = "ppt") %>% 
  select(ppt, RDE, FE)

RDE_FE1 <- RDE_FE1 %>% 
  mutate(difference = FE - RDE)

FEvsRDE <- t.test(RDE_FE1$FE, RDE_FE1$RDE, paired = TRUE)
tidy(FEvsRDE)
```

```{r RDE vs FE analyses non-para, message = FALSE}
# Wilcoxon Signed rank test for RDE vs FE
FE_RDEtest <- wilcox.test(RDE_FE1$FE, RDE_FE1$RDE, paired = TRUE)
tidy(FE_RDEtest)
```

```{r RDE vs FE2, message = FALSE, echo=FALSE, results='hide'}
mean(RDE_FE1$difference)
sd(RDE_FE1$difference)

mean(RDE_FE1$RDE)
sd(RDE_FE1$RDE)
mean(RDE_FE1$FE)
sd(RDE_FE1$FE)
```

In fact, not only were familiarity effects more pronounced than reverse distance effects at the group level (t(55) = 8.83, p < .001) [V = 1555, p < .001], but also at the individual level, 94.64% of participants (53 out 56) displayed a familiarity effect larger than their reverse distance effect.

This effect remained even when only considering those who displayed a RDE:

```{r only those who display a RDE, message = FALSE, echo=TRUE}
RDEgreaterthan0only <- RDE_FE1 %>% 
  filter(RDE > 0)

```

```{r show these results, message = FALSE, echo=FALSE}
### with only those who show a RDE
FEvsRDE2 <- t.test(RDEgreaterthan0only$FE, RDEgreaterthan0only$RDE, paired = TRUE)
tidy(FEvsRDE2)

wilcox.test(RDEgreaterthan0only$FE, RDEgreaterthan0only$RDE, paired = TRUE)

```



Notably, 28.57% of participants (16 out 56) had a difference smaller than 0 ms, indicating they did not display a reverse distance effect but rather processed non-consecutive sequences faster than consecutive sequences (i.e., a standard distance effect).

## Arithmetic fluency analysis

```{r load in arithmetic data again, include=FALSE, echo=FALSE, warning=FALSE}
arithmetic_data <- read_excel("data.xlsx", sheet = "TTR")

TTR <- arithmetic_data %>% 
  select(ppt, total)

seqdetails2 <- type %>% 
  mutate(type2 = case_when(direction == "A" & distance == "C" ~ "AC",
                           direction == "A" & distance == "N" ~ "ANC",
                           direction == "D" & distance == "C" ~ "DC",
                           direction == "D" & distance == "N" ~ "DNC")) %>% 
  select(seq, type2)

arith_df <- order_df %>% 
  left_join(seqdetails2, by = "seq") %>% 
  select(ppt, seq, RT, type = type2) %>% 
  group_by(ppt, type) %>% 
  summarise(medRT = median(RT)) %>% 
  ungroup() %>% 
  pivot_wider(names_from = type, values_from = medRT) %>% 
  left_join(FE_df, by = "ppt") %>% 
  left_join(TTR, by = "ppt")

all <- order_df %>% 
  group_by(ppt) %>% 
  summarise(all = median(RT)) %>% 
  ungroup()

arith_df <- arith_df %>% 
  left_join(all, by = "ppt")
```

Finally, we explored the association between order verification and arithmetic fluency. Correlations were calculated between overall arithmetic scores and response times for six sequence types: ascending consecutive, ascending non-consecutive, descending consecutive, descending non-consecutive, familiar, and unfamiliar:

```{r new corr matrix, message = FALSE, results='hide', echo=FALSE, warning=FALSE}
# run new correlation matrix
#arith_df <- arith_df %>% 
#  select(-ppt)

arith_df <- arith_df %>% 
  rename(TTR = total, familiar = high, unfamiliar = low) %>% 
  select(-FE)

arith_df <- arith_df %>% 
  select(-ppt)

arith_df <- arith_df %>% 
  rename("Unfamiliar" = unfamiliar, "Familiar" = familiar, "Arithmetic" = TTR)

cor <- round(cor(arith_df), 2)
p.mat <- cor_pmat(arith_df)

```

```{r new corr matrix2, message = FALSE, echo=FALSE, results='hide', warning=FALSE}
# calculate p values
tidy(p.mat)


# Create the ggcorrplot
my_plot <- ggcorrplot(cor, method = "square", hc.order = TRUE, type = "lower", 
                      title = "", legend.title = "Pearson\n corr",
                      lab = TRUE, lab_col = "black", lab_size = 3,
                      ggtheme = theme_classic, outline.color = "white",
                      colors = c("deepskyblue1", "white", "#FF6A6A"), p.mat = p.mat, insig = "blank")

# Create a text box with the caption
my_caption <- textGrob("*All displayed correlations are statistically significant (p < .05)", 
                       gp = gpar(fontsize = 10), 
                       x = 0.5, y = 0.95)
```

```{r show corr matrix, message = FALSE, echo=FALSE, results='hide', out.width="95%"}
grid.arrange(my_plot, my_caption, nrow = 2, heights = c(4,0.5))
```

Overall, there was a moderate negative correlation between response times and arithmetic scores, suggesting faster order verification is associated with higher arithmetic fluency. 

This association was strongest for descending consecutive sequences and weakest for ascending non-consecutive sequences. However, there appeared no obvious pattern of associations between specific sequence types and arithmetic scores, with correlations ranging between -0.41 and -0.53. 

Notably, response times for familiar (-0.46) and unfamiliar (-0.47) sequences showed similar associations with arithmetic fluency.