---
title: "Attention is spotaneously biased towards statistical regularities." 
shorttitle: ""
 
author:
  - name            : "Alexandra Vlassova"
    affiliation     : "1"
    corresponding   : yes
    address         : .
    email           : "a.vlassova@donders.ru.nl"

  - name            : "Eelke Spaak"
    affiliation     : "1"

  - name            : "Floris P. de Lange"
    affiliation     : "1"

affiliation:
  - id              : "1"
    institution     : "Donders Centre for Cognitive Neuroimaging, Radboud University Nijmegen"
  
#abstract: |
#  Will eventually have words that go here. 
   
# keywords          : "keywords"
# wordcount         : "X" 
bibliography      : ["references.bib"]

floatsintext      : no
figurelist        : no
tablelist         : no
footnotelist      : no
linenumbers       : no
mask              : no
draft             : no

documentclass     : "apa6"
classoption       : "man"
output: papaja::apa6_word  
csl: royal-society-open-science.csl
  
---

```{=html}
<style>
p.caption {
  font-size: 0.6em;
}
</style>
```

```{r setup, include = FALSE}
# Load packages ----------------------------------------------------------------
library(knitr)
library(papaja)
library(tidyverse)
library(rstatix)
library(BayesFactor)
library(e1071)
library(ggpubr)
library(patchwork)
library(cowplot)
library(lattice)
library(scales)
library(png)
library(afex)
library(emmeans)


# Compile settings -------------------------------------------------------------
opts_chunk$set( 
	fig.path = "figures/",
	message = FALSE,
	cache = FALSE,
	cache.extra = knitr::rand_seed,
	dev = "png",
	dpi=300 
	) 

session_info <- sessionInfo()  
set.seed(64)

# Rounding functions -----------------------------------------------------------
dc2 <- function(x, d = 2) sprintf(paste0("%1.", d, "f"), x)
dc1 <- function(x, d = 1) sprintf(paste0("%1.", d, "f"), x)
dc0 <- function(x, d = 0) sprintf(paste0("%1.", d, "f"), x)

r_refs("references.bib")
```

# Summary
In this pre-registered replication study, we reassessed the spontaneous bias of attention towards task-irrelevant statistical regularities reported by Zhao, Al-Aidroos, and Turk-Browne (Psychological Science, 24(5), 667–677, 2013) using a larger sample size. We additionally extended the study to investigate the potential role of learning, doubling the duration of the experiment to determine whether the attentional bias towards regularities is reduced with continued exposure. While we could not directly replicate Zhao et al.’s central findings, our results provide a conceptual replication: we found that performance on a visual search task was better in a location containing task-irrelevant statistical regularities, supporting their conclusion that attention is spontaneously biased towards regularities in the environment. We further found that the attentional bias effect diminished with continued exposure to the regularities. We conclude that the spontaneous bias of attention towards task-irrelevant regularities may be transient, and propose that learning of the regularities may play a key role in moderating attentional selection.  

# Introduction

Our environment is complex and filled with a vast amount of information, yet we are limited in our capacity to process it all. We can optimise how we deploy our limited cognitive resources by effectively allocating attention to prioritise the extraction of relevant information, while ignoring distracting information. Much of the incoming information is highly structured, containing regularities in objects and events that occur across space and time. Statistical learning refers to the ability to extract these regularities from the world around us, which in turn allows us to more effectively select and process information [@Schapiro2015; @Frost2019]. Despite a wealth of studies across sensory and cognitive domains demonstrating that statistical regularities can be learned automatically [@Fiser2001; @Saffran1996; @Saffran1999; @Kim2009; @Perruchet2006; @Batterink2019; @Spaak2020] – even when they are entirely irrelevant to the task at hand [@Turk-Browne2005; @Duncan2020; @Musz2015] – the relationship between attention and statistical learning remains poorly understood.

Attentional selection is typically considered to be directed by the interaction between goal-driven top-down attention and stimulus-driven bottom-up attention [@Theeuwes2010; @Itti2001]. However, this dichotomy fails to account for cases in which biases in attention cannot be explained by current goals or stimulus salience alone. An example of this is the study by Zhao et al. [@Zhao2013], who proposed that task-irrelevant regularities might automatically attract attention, allowing for statistical learning to occur more rapidly. They presented sequences of abstract shapes in four fixed locations, with the order of the sequences statistically structured in one of the locations and random in the others. These shape sequences were occasionally interrupted at random intervals by a visual search array. In a visual search task, participants search for a predeﬁned target stimulus among distractors and are asked to respond as fast as possible. The purpose of the search task in Zhao et al. [@Zhao2013] was to assess attentional allocation across the four sequence locations. If, as the authors hypothesised, spatial attention is biased towards statistical regularities, then discrimination on the visual search task should be faster when the visual search targets appear in the structured location than in the random location, as attention would be more likely to have already been directed at that location. They found that response times were indeed faster to targets that were presented in the location containing the structured sequence, despite the sequences not having any predictive information as to the target location and there being no benefit of attending to the structured stream for finding the target on the visual search task. Participants were also unable to correctly identify the location of the structured stream, with most reporting that they had not noticed any regularities during the task, indicating that the attentional bias effect could not be simply explained by intentional orienting to the structured location. Zhao et al. therefore concluded that attention was spontaneously biased towards spatial locations containing task-irrelevant statistical regularities, as evidenced by improved performance on the visual search task at that location. They proposed that the attentional bias towards regularities may be due to prioritising attentional resources on stable and learnable aspects of the environment.   

While Zhao et al.’s findings offer an impactful insight into the effect of statistical learning on attentional allocation, their sample size and effect size were relatively modest, rendering the results somewhat tentative. To address these limitations, we performed an exact replication of Zhao et al.’s Experiment 1 [@Zhao2013] in a sample of 91 participants. We additionally doubled the duration of the experiment to increase exposure to the regularities in order to investigate whether learning plays a moderating role in the attentional bias effect. If the attentional bias towards regularities serves to increase the automatic learning of structure, one would hypothesize the effect to be transient, only present while the regularities are learnt. If, on the other hand, humans show a more general bias towards structure, independent of the learning process, one would expect the effect to be stable. 
 
# Methods 

```{r import_data, include = FALSE}
setwd("/Users/alya/Documents/PROJECTS/PREDATT/EXPERIMENT1/ANALYSIS/")

# Import task data -------------------------------------------------------------
data_import <- read_csv(
  file = "./DATA/PA_EXP1_REP.csv",
  col_types = cols("f", "i", "f", "f", "f", "i", "f", "f", "f", "d")
  ) %>%
  mutate(correct = ifelse(search_rotation == response, 1, 0)) %>%
  mutate(target_location = factor(
    ifelse(search_location == structured_location, "structured", "random"),
    levels = c("structured", "random")
  ))

# Import castor data -----------------------------------------------------------
castor_import <- read_csv(
  file = "./DATA/castor/main.csv",
  col_types = cols_only("i", "f", "d", "f", "c", "i", "i", "f", "i", "f", "f")
  ) %>%
  rename(
    id = "Record Id",
    subject = info_general_sonaID_1,
    dob = info_general_year,
    gender = info_general_gender
  ) %>%
  filter(!id %in% c(110015, 110016, 110017, 110018, 110019, 110048)) %>%
  select(subject, dob, gender, q_regularity, q_regularityloc) %>% # remove technical error records
  transform(
    q_regularityloc =
      ifelse(q_regularityloc == "Left", 1,
        ifelse(q_regularityloc == "Right", 2,
          ifelse(q_regularityloc == "Top", 3,
            ifelse(q_regularityloc == "Bottom", 4, NA)
          )
        )
      )
  )

# Import familiarity test data -------------------------------------------------
mem_data <- read_csv(file = "./DATA/PA_EXP1_MEM.csv", col_types = cols("f", "i", "f", "f", "d"))
```

## Sample size
We initially preregistered this experiment (including hypotheses, experimental designs, and analysis plans) prior to data collection on OSF (<https://osf.io/ngv5s>). This article received results-blind in-principle acceptance (IPA) at Royal Society Open Science. Following IPA, the accepted Stage 1 version of the manuscript, not including results and discussion, was also preregistered on the OSF (<https://osf.io/zkthf>). This preregistration was performed after data analysis.

We calculated that we would need to collect data from 64 participants to ensure \>95% power to find an effect equal to or larger than Zhao et. al's [@Zhao2013] reported effect size. However, as we also intended to look at the correlation between accuracy on a surprise Familiarity test and the differences in RT on the visual search task, we determined that we would require 84 subjects to ensure a sufficient sample size for a Pearson correlation with an alpha of .05, 80% power, and a medium effect size (*r* = .3) using a two-tailed analysis. These sample sizes were calculated using G\*Power [Version 3.1.9.6; @Faul2009] and pre-registered.

## Participants

```{r demographics, include = FALSE}
age <- castor_import %>%
  mutate(age = 2020 - dob) %>% 
  summarise(n = n(), mean_age = mean(age), sd = sd(age), min = min(age), max = max(age), .groups = "drop")

gender <- castor_import %>%
  group_by(gender) %>%
  summarise(count = n_distinct(subject), .groups = "drop")
```
 
A total of 91 subjects were recruited for this experiment. Data from 5 participants were removed due to a technical issue which affected the precision of the recorded response times. One additional participant was excluded as they were unable to finish the experiment within the allocated time window. These participants were replaced in order to attain the pre-registered goal of N = 84 participants. In total, we collected `r age$n` complete datasets (`r gender$count[1]` male; `r age$min` - `r age$max` years of age, *M* = `r dc0(age$mean_age)`, SD = `r dc1(age$sd)`).

Participants were recruited through the Radboud University participant pool and had normal or corrected-to-normal vision. They provided written informed consent and were compensated for their time at a rate of 8 Euros per hour, regardless of performance. The study was approved by the local ethics committee (CMO Arnhem -- Nijmegen, Radboud University Medical Center) under the general ethics approval ("Imaging Human Cognition", CMO 2014/288).

## Apparatus

Participants were seated on a height-adjustable chair at a distance of approximately 60 cm from a 24" BenQ XL2420T monitor, with a resolution of 1920 x 1090 and a 60 Hz refresh rate. Stimuli were presented using Psychtoolbox Version 3 [@Brainard1997] for MATLAB 2016a (The MathWorks, Natick, MA) on a Dell Precision Tower 5810 desktop computer running Windows 7 OS.

## Stimuli

Zhao et al. kindly shared their original shape stimuli set with us, and unless otherwise noted, all other details were matched as closely as possible to the specifications outlined in Zhao et al. [@Zhao2013], Experiment 1.

Shapes were presented on a white background in four locations: above, below, left and right of the central fixation point (Figure 1A). Each shape subtended 3.3^o^ and was centred 5.1^o^ from the fixation point, with each location outlined with a black square outline which remained on the screen. The stimuli shown at each location were selected from four separate sets of shapes each consisting of nine shapes, and were presented simultaneously on each trial. Each set was presented 100 times (50 times in Block 1 as in Zhao et al., and an additional 50 times in Block 2, as per our pre-registered extension to the design), resulting in 900 total trials (9 shapes per set x 100 repetitions).

One of the four stimulus sets was presented in the “structured” location, with both the chosen set and the location counterbalanced across participants. In the structured location, the nine shapes were randomly grouped for each participant into three triplets, so that while the order in which the triplets were presented was shuffled (with the constraint that no triplet could be repeated twice in a row), the order of the shapes within each triplet remained constant throughout the experiment. In the three “random” locations, the order in which the shapes from each set were presented were randomised, with the constraint that no shape could be repeated consecutively.

The sequence of shape presentations was occasionally interrupted by a visual search array (Figure 1B). Each array contained one target “T” shape, and three “L”-shaped distractors, each rotated so that the long line was horizontal and the short line would be pointing leftwards or rightwards. The pointing directions were counterbalanced so that one distractor pointed in the same direction as the target and the other two distractors pointed in the opposite direction. The distractor’s horizontal lines were offset by 20% from either the top or bottom of the vertical line. The visual search shapes were presented centered in the black square outline that remained on the screen throughout both the shape and visual search presentation. The visual search targets appeared with equal frequency at each of the of the four locations, and the sequences of target locations and directions were randomly shuffled within each block so that neither could be predicted from previous trials. The target location was therefore uninformative as to which location contained the structured stimulus set, and attending to the structured location would not benefit finding the target on the visual search task.

## Procedure

As in Zhao et al. [@Zhao2013], participants were asked to fixate while attending to the four stimulus locations, which were marked with a black square outline. On each shape trial, one shape from each of the four corresponding stimulus sets was presented in each square for 750 ms, followed by an interstimulus interval (ISI) of 750 ms, during which only the square outline and fixation point remained on the screen. There were 900 such shape trials, with 160 search trials randomly interspersed with the restriction that a search trial could not occur twice in a row, resulting in an average of 6.6 shape trials (*SD* = 4.8) preceding each search trial. On average, participants saw 6.6 shape trials (*SD* = 0.04; range = 6.4 – 6.6 trials) before each search trial, with search trials separated by a minimum of 1.7 (*SD* = 0.5) and a maximum of 28.7 (*SD* = 4.8) shape trials.

The visual search arrays were presented for 750 ms, followed by an ISI of at least 750 ms, during which only the square outline and fixation point remained on the screen. Participants were instructed to responded as quickly and accurately as possible, using a standard keyboard to indicate whether the target pointed left (“1” key) or right (“0” key). If the participant didn’t respond during the array presentation or ISI, the ISI was extended until a response was made. After completing Block 1 (450 shape trials and 80 visual search trials), participants were given the opportunity to take a short break prior to completing the second half of the trials (Block 2). Prior to beginning the experiment, participants were shown 5 random shape trials followed by 1 search trial, so that they could practice the task and clarify the instructions.

Once participants had completed both blocks, their learning of the regularities was probed with a surprise two-alternative forced-choice test, in which they were asked to distinguish between shape sequences that differed solely in the familiarity of their higher-order temporal structure. On each trial, participants were shown two sequences of three shapes presented at fixation, with a 1000 ms pause between sequences. In one sequence, shapes from a triplet that had appeared repeatedly in the structured stream during the main task, were presented one at a time for 750 ms with a 750 ms ISI. The other sequence was a “foil,” composed by selecting one shape from each of the three triplets such that the position of the shape within its original triplet was preserved (i.e. if the triplets in the structured location were shapes “ABC,” “DEF,” and “GHI,” foils would be “AEI,” “DHC,” “GBF”). Participants used the “1” and “2” keys to indicate which sequence seemed more familiar to them based on what they experienced in the main task. Each triplet was tested against each foil twice, resulting in a total of 18 trials. Trial order was randomised and sequence presentation order counterbalanced across trials.

After the test, participants were asked to complete a brief questionnaire, where they were asked what they thought the experiment was testing, whether they followed a strategy of attending to a specific location (and if so, what that strategy was), whether they noticed any regularities in the order the shapes were presented in, and how confident they were in their performance on the familiarity test. Participants were then debriefed and told about the triplets in the main task, after which they were asked to choose the location they thought contained the structured stream, and rate their confidence in this response.

We followed the details of the procedure used by Zhao et al. (Experiment 1) as closely as possible, with Block 1 being a direct replication. The key difference between our designs was that we doubled the number of trials in order to look at potential changes in the attentional bias effect with learning.

## Data Analysis

All data and code (including the processed data and data-analysis files in the programming languages R, and experimental scripts in MATLAB) are openly available on the Donders Institute for Brain, Cognition and Behavior repository (<https://data.donders.ru.nl/collections/di/dccn/DSC_3018028.05_239>).

To ensure high data quality, our pre-registered exclusion criteria were to reject participants whose overall accuracy fell more than two standard deviations below the group mean. As per Zhao et al. and our pre-registered criteria, we included only correct responses in the RT analysis, and removed trials with RTs greater than 3 standard deviations from the mean in each condition.  

In order to be able to directly compare our results with those of Zhao et al. [@Zhao2013], differences in performance (RTs and error rates) on the search task between the structured and random locations were first examined using paired t-tests. As per our pre-registered plans, we conducted a paired t-test on the differences in RT between blocks to determine whether the potential RT facilitation in Block 1 is reduced during Block 2, as well as on the relative RT differences (\% change from Block 1 to Block 2), to ensure that potential differences do not simply reflect a training-driven proportional reduction in RTs. We additionally conducted a 2 (structured or random locations) x 2 (block) repeated-measures analysis of variance (rmANOVA), on both RTs and accuracy to investigate potential interactions between the attentional bias effect and block.  

In order to determine whether a preference towards one of the locations improved performance at that location independently from the presence of regularities, we examined performance on the search task based on participants’ reports of where they believed (or guessed) the location of the structured stream was. That is, rather than comparisons between the *actual* structured and random locations, we instead used an rmANOVA to compare performance between the location participants *thought* contained the structured stream with the others. We additionally used the accuracy of the reported structured location in the debrief as a between-subjects factor in an rmANOVA to check whether there were differences in the performance facilitation effects between subjects who accurately reported the location of the structured stream and those who did not. Individual contrasts were tested using Bonferroni-adjusted pairwise t-tests.

As in Zhao et al. [@Zhao2013], we used a one sample t-test to determine whether accuracy on the 2AFC familiarity test was above chance levels. Following our pre-registered plans, we additionally investigated whether there was an association between learning of the regularities (as measured by the 2-AFC test) and the attentional capture effect in each block, as well as the difference in this effect between blocks using Pearson correlation tests. We conducted an equality of proportions test to determine whether participants who reported being aware of the regularities in the debrief were also more likely to correctly identify the location containing the structured stream, and we performed rmANOVAs to determine if there were differences in performance on the familiarity test between those who reported being aware of the regularities and those who did not, as well as between those who correctly reported the location of the regularities and those who did not.


In addition to looking at differences in RTs and accuracy in isolation, to address potential differences in speed-accuracy tradeoffs we chose to calculate linear integrated speed accuracy scores [LISAS; @Vandierendonck2017], which allowed us to combine RT and accuracy into a single score. This measure obtains similar weights of the two components by weighting the accuracy with the ratio of the RT and accuracy standard deviations as follows:

$$LISAS_j\ =\ RT_j\ +\ \left(\frac{S_{RT\ }}{S_{PE}}\right)\ \times PE_j$$

where RT~*j*~ is the participant's mean RT in condition *j*, PE~*j*~ is the participant's proportion of errors in condition *j*, S~*RT*~ is the participant's overall RT standard deviation, and S~*PE*~ is the participant's overall PE standard deviation. We repeated the analyses described above using the LISAS scores.

All tests were two-tailed since the original study had a small sample and effect sizes, and we thought it possible that the effect could go in the opposite direction in our replication. We report p-values (\(\alpha\) = 0.05) and Bayes Factors (BF~10~) for all analyses. All analyses were performed using R (v. 4.0.2).


# Results 

```{r clean_accuracy, message=TRUE, warning=TRUE, include=FALSE}
below_chance <- data_import %>%
  group_by(subject) %>%
  summarise(accuracy = mean(correct) * 100, .groups = "drop") %>%
  filter(accuracy < 50)

# remove subjects performing below chance overall & extreme RTs
data_raw <- filter(data_import, !subject %in% below_chance$subject) %>%
  filter(RT < 10)

# remove subjects with accuracy >2SD below the group mean
acc_cutoff <- data_raw %>%
  group_by(subject) %>%
  summarise(accuracy = mean(correct) * 100, .groups = "drop") %>%
  mutate(cutoff = mean(accuracy) - 2 * sd(accuracy)) %>%
  filter(accuracy < cutoff)

data_raw <- data_raw %>% mutate(exc = ifelse(subject %in% acc_cutoff$subject, 1, 0))
```

Three participants had accuracy rates well below chance levels on the visual search task (`r dc0(below_chance[1,2])`% - `r dc0(below_chance[3,2])`% correct; chance level = 50%), which we consider to reflect a failure to follow task instructions rather than genuine task performance. To avoid artificially lowering the cutoff point, we first excluded these subjects before calculating the accuracy cutoff according to our pre-registered plan (accuracy \>2SD below the group mean). Based on this criterion, we excluded an additional `r nrow(acc_cutoff)` participants who had overall accuracy levels below the cutoff of `r dc0(acc_cutoff$cutoff[1])`%, leaving a total of `r n_distinct(filter(data_raw, exc==0)$subject)` participants.


```{r skew_sample, message=TRUE, warning=TRUE, include=FALSE}
sample_size <- filter(data_raw, exc==0) %>%
  group_by(subject) %>%
  summarise(accuracy = mean(correct), .groups = "drop") %>%
  summarise(
    skew = skewness(accuracy),
    n = -5.19 * skew^4 - 1.44 * skew^3 + 5.65 * skew^2 + 14.75 * skew + 70.07
  )
 
```

While the distributions of both the accuracy and the differences between conditions were neither normal nor symmetrical, our sample size was large enough to obtain stable estimates (given the skewness level, we would require a minimum sample size of `r round(sample_size$n)` [@Piovesana2018]), allowing us to follow the pre-registered plans to use the parametric tests as originally reported in Turk-Browne.

```{r clean_RT, message=TRUE, warning=TRUE, include=FALSE}
# calculate RT cutoffs (>3SD from means, correct trials only)
rt_cutoff <- filter(data_raw, correct==1) %>% 
  group_by(subject, block, target_location) %>%
  summarise(mRT = mean(RT), sdRT = sd(RT), .groups = "rowwise") %>%
  mutate(max = mean(mRT) + (3 * sdRT)) 
 
# flag outliers in main df
data_raw <- left_join(data_raw, rt_cutoff) %>%
  mutate(flag = ifelse(RT > max, 1, 0))

# summarise outliers
rt_outliers <- filter(data_raw, exc==0 & correct==1)  %>%
  group_by(block, target_location, subject) %>%
  summarise(out = sum(flag), n = n(), p = out / n * 100, .groups = "drop_last") %>%
  summarise(percent = mean(p), sd = sd(p), max = max(p), .groups = "drop")
```
 
As per Zhao et al. and our pre-registered criteria, we included only correct responses in the RT analysis, and removed trials with RTs greater than 3 standard deviations from the mean in each condition. An average of `r dc1(rt_outliers$percent[1])`% (*SD* = `r dc1(rt_outliers$sd[1])`%) of trials in the structured condition and `r dc1(rt_outliers$percent[2])`% (*SD* = `r dc1(rt_outliers$sd[2])`%) in the random condition were identified as outliers due to an RT slower than the cutoff and were removed for subsequent RT analyses.

```{r means, message=TRUE, warning=TRUE, include=FALSE}
accuracy_means <- data_raw %>%
  group_by(subject, target_location, block, exc) %>%
  summarise(accuracy = mean(correct)*100, .groups = "drop") 

data_means <- filter(data_raw, flag==0 & correct==1) %>%
  group_by(subject, target_location, block, exc) %>%
  summarise(rt = mean(RT), .groups = "drop") %>% 
  left_join(accuracy_means) 
```

## Direct Replication

```{r block1, message=TRUE, warning=TRUE, include=FALSE}
rand <- filter(data_means, block==1 & exc==0 & target_location=="random")
struct <- filter(data_means, block==1 & exc==0 & target_location=="structured")

# RT ---------------------------------------------------------------------------
b1_rt_TT <- t.test(rt ~ target_location, filter(data_means, block==1 & exc==0), paired = TRUE)
b1_rt_CD <- cohens_d(filter(data_means, block==1 & exc==0), rt ~ target_location, paired = TRUE)
b1_rt_BF <- ttestBF(rand$rt, struct$rt, paired = TRUE)

# Accuracy ---------------------------------------------------------------------
b1_accuracy_TT <- t.test(accuracy ~ target_location, filter(data_means, block==1 & exc==0), paired = TRUE)
b1_accuracy_CD <- cohens_d(filter(data_means, block==1 & exc==0), accuracy ~ target_location, paired = TRUE)
b1_accuracy_BF <- ttestBF(rand$accuracy, struct$accuracy, paired = TRUE)
```

Zhao et al. found that RTs were significantly faster when the visual search targets appeared in the structured location, and did not report any significant differences in accuracy [@Zhao2013]. We first looked at performance in Block 1, as this was a direct replication of Zhao et al.'s Experiment 1. Unlike Zhao et al., paired t-tests revealed no significant differences in RTs between trials in which the visual search target appeared in the structured location and trials where the target appeared in one of the three random locations (`r apa_print(b1_rt_TT)$full_result`, $d$ = `r b1_rt_CD$effsize`, *BF~10~* = `r as.vector(b1_rt_BF)`; Figure 2A). However, we found that subjects had significantly higher accuracy on the visual search task when the target appeared in the structured location compared to when it was presented in the random location (`r apa_print(b1_accuracy_TT)$full_result`, $d$ = `r b1_accuracy_CD$effsize`, *BF~10~* = `r as.vector(b1_accuracy_BF)`; Figure 2B).

```{r block1_noexc, message=TRUE, warning=TRUE, include=FALSE}
rand <- filter(data_means, block==1 & target_location=="random")
struct <- filter(data_means, block==1 & target_location=="structured")
 
# RT ---------------------------------------------------------------------------
b1_rt_TT_noexc <- t.test(rt ~ target_location, filter(data_means, block==1), paired = TRUE)
b1_rt_CD_noexc <- cohens_d(filter(data_means, block==1), rt ~ target_location, paired = TRUE)
b1_rt_BF_noexc <- ttestBF(rand$rt, struct$rt, paired = TRUE)

# Accuracy ---------------------------------------------------------------------
b1_accuracy_TT_noexc <- t.test(accuracy ~ target_location, filter(data_means, block==1), paired = TRUE)
b1_accuracy_CD_noexc <- cohens_d(filter(data_means, block==1), accuracy ~ target_location, paired = TRUE)
b1_accuracy_BF_noexc <- ttestBF(rand$accuracy, struct$accuracy, paired = TRUE)
```

Since Zhao et al. did not report excluding any participants based on their accuracy, we repeated the analyses without excluding participants on the basis of the 2SD criterion, and found that this did not change the pattern of significance (RT: `r apa_print(b1_rt_TT_noexc)$full_result`, $d$ = `r b1_rt_CD_noexc$effsize`, *BF~10~* = `r as.vector(b1_rt_BF_noexc)`; accuracy: `r apa_print(b1_accuracy_TT_noexc)$full_result`, $d$ = `r b1_accuracy_CD_noexc$effsize`, *BF~10~* = `r as.vector(b1_accuracy_BF_noexc)`).

In summary, while we did not directly replicate the behavioural pattern reported by Zhao et al., we did find that performance on the visual search task was improved when the targets appeared in the structured location, supporting Zhao et al.'s conclusion that attention is spontaneously biased to locations that contain statistical regularities, even when these regularities are task irrelevant. However, while Zhao et al. found that performance was improved by speeding up response times, we found that the performance improvements were particularly visible in the accuracy rates. A potential reason for this difference may be as simple as subtle differences in task instructions, which have been previously shown to modulate speed-accuracy trade-offs [@vanVeen2008].

```{r lisasb1, message=TRUE, warning=TRUE, include=FALSE}
lisas <- filter(data_raw, flag==0 & exc==0) %>%
  group_by(subject) %>%
  mutate(S_RT = sd(RT)) %>%
  group_by(subject, target_location, block) %>%
  summarise(
    RTj = mean(RT),
    PEj = 1 - mean(correct),
    S_RT = max(S_RT),
    .groups = "drop"
  ) %>%
  group_by(subject) %>%
  mutate(S_PE = sd(PEj)) %>%
  group_by(subject, target_location, block) %>%
  summarise(score = RTj + (S_RT / S_PE) * PEj, .groups = "drop")
  
rand <- filter(lisas, block==1 & target_location=="random")
struct <- filter(lisas, block==1 & target_location=="structured")

b1_lisas_TT <- t.test(score ~ target_location, filter(lisas, block==1), paired = TRUE)
b1_lisas_CD <- cohens_d(filter(lisas, block==1), score ~ target_location, paired = TRUE)
b1_lisas_BF <- ttestBF(rand$score, struct$score, paired = TRUE)  
```

We repeated the same paired t-test analysis using the LISAS, and found that scores were significantly smaller in the structured location than in the random locations in the first block (`r apa_print(b1_lisas_TT)$full_result`, $d$ = `r b1_lisas_CD$effsize`). A Bayes factor analysis provided strong evidence (*BF~10~* = `r as.vector(b1_lisas_BF)`) to support our initial conclusion that performance on the visual search task was better when the visual search targets appeared in the structured location.


## Attentional bias towards regularities decreases with exposure.

We next considered whether the magnitude of the attentional bias effect is constant throughout the experiment, or whether it decreases with exposure to the regularities. 

```{r block2, message=TRUE, warning=TRUE, include=FALSE}
rand <- filter(data_means, block==2 & exc==0 & target_location=="random")
struct <- filter(data_means, block==2 & exc==0 & target_location=="structured")

# RT ---------------------------------------------------------------------------
b2_rt_TT <- t.test(rt ~ target_location, filter(data_means, block==2 & exc==0), paired = TRUE)
b2_rt_CD <- cohens_d(filter(data_means, block==2 & exc==0), rt ~ target_location, paired = TRUE)
b2_rt_BF <- ttestBF(rand$rt, struct$rt, paired = TRUE)

# Accuracy ---------------------------------------------------------------------
b2_accuracy_TT <- t.test(accuracy ~ target_location, filter(data_means, block==2 & exc==0), paired = TRUE)
b2_accuracy_CD <- cohens_d(filter(data_means, block==2 & exc==0), accuracy ~ target_location, paired = TRUE)
b2_accuracy_BF <- ttestBF(rand$accuracy, struct$accuracy, paired = TRUE)
```

When inspecting visual search behaviour in the second part of the experiment (Block 2), we found that there were no significant differences in visual search performance when targets appeared at a structured compared to random location, neither in RT (`r apa_print(b2_rt_TT)$full_result`, $d$ = `r b2_rt_CD$effsize`, *BF~10~* = `r as.vector(b2_rt_BF)`), nor in accuracy (`r apa_print(b2_accuracy_TT)$full_result`, $d$ = `r b2_accuracy_CD$effsize`, *BF~10~* = `r as.vector(b2_accuracy_BF)`; Figure 2B).

```{r rt_diff, message=TRUE, warning=TRUE, include=FALSE}
diffs <- filter(data_means, exc==0) %>%  
  pivot_wider(names_from = target_location, values_from = c(rt, accuracy)) %>% 
  mutate(rt = rt_structured - rt_random, 
         accuracy = accuracy_structured - accuracy_random)

diff_means <- diffs %>% 
  group_by(block) %>% 
  summarise_at(c("accuracy","rt"), list(m = mean, sd = sd))

diff_rt_TT <- t.test(rt ~ block, diffs, paired = TRUE)
diff_rt_CD <- cohens_d(diffs, rt ~ block, paired = TRUE)
diff_rt_BF <- ttestBF(filter(diffs, block==1)$rt, filter(diffs, block==2)$rt, paired = TRUE)
```

To directly compare whether the effect of target location on performance was different between the two blocks, we conducted a paired samples t-test on the mean differences in RT and accuracy between conditions in Block 1 compared to the mean differences in Block 2 (i.e. is there an effect of block on the difference between conditions?). On average, participants had a significantly larger RT difference between the two conditions in Block 2  (*M* = `r diff_means$rt_m[2]` s difference between conditions (structured - random), *SD* = `r diff_means$rt_sd[2]`) compared to Block 1 (*M* = `r diff_means$rt_m[1]` s, *SD* = `r diff_means$rt_sd[1]`; `r apa_print(diff_rt_TT)$full_result`, $d$ = `r diff_rt_CD$effsize`), though the evidence for this was only anecdotal (*BF~10~* = `r as.vector(diff_rt_BF)`). 

```{r rt_reldiff, message=TRUE, warning=TRUE, include=FALSE}
rt_rel <- filter(data_means, exc==0) %>%  
  pivot_wider(id_cols = c(subject, target_location), names_from = block, values_from = rt) %>% 
  mutate(diff = `1` - `2`, reldiff = diff/`1`*100) 

rt_rel_mean <- rt_rel %>% 
  group_by(target_location) %>% 
  summarise(m = mean(reldiff), sd = sd(reldiff), .groups = "drop")

rand = filter(rt_rel, target_location == "random")
struct = filter(rt_rel, target_location == "structured")

rel_rt_TT <- t.test(reldiff ~ target_location, rt_rel, paired = TRUE) 
rel_rt_CD <- cohens_d(rt_rel, reldiff ~ target_location, paired = TRUE)
rel_rt_BF <- ttestBF(struct$reldiff, rand$reldiff, paired = TRUE)
```

To determine whether this RT facilitation effect could be explained by a training-driven proportional reduction in RTs, we examined the relative changes in RT differences between blocks. On average, participants were `r dc1(rt_rel_mean$m[1])`% (*SD* = `r dc1(rt_rel_mean$sd[1])`%) faster in the structured condition in Block 2 compared to their performance in the structured condition in Block 1, and `r dc1(rt_rel_mean$m[2])`% (*SD* = `r dc1(rt_rel_mean$sd[2])`%) faster in the random condition. The improvement in performance between blocks was significantly larger in the structured than in the random condition (`r apa_print(rel_rt_TT)$full_result`, $d$ = `r rel_rt_CD$effsize`, *BF~10~* = `r as.vector(rel_rt_BF)`), indicating that RTs improved in the structured condition beyond what can be accounted for by a general training effect or baseline differences between conditions.

```{r accuracy_diff, message=TRUE, warning=TRUE, include=FALSE}
diff_accuracy_TT <- t.test(accuracy ~ block, diffs, paired = TRUE)
diff_accuracy_CD <- cohens_d(diffs, accuracy ~ block, paired = TRUE)
diff_accuracy_BF <- ttestBF(subset(diffs, block==1)$accuracy, subset(diffs, block==2)$accuracy, paired = TRUE)
```

For accuracy, the mean difference between conditions in Block 1 (*M* = `r diff_means$accuracy_m[1]` % difference between conditions (structured - random), *SD* = `r diff_means$accuracy_sd[1]`) was slightly reduced in Block 2 (*M* = `r diff_means$accuracy_m[2]` %, *SD* = `r diff_means$accuracy_sd[2]`), but the reduction of this difference between blocks was not statistically significant (`r apa_print(diff_accuracy_TT)$full_result`, $d$ = `r diff_accuracy_CD$effsize`, *BF~10~* = `r as.vector(diff_accuracy_BF)`). 

```{r lisas, message=TRUE, warning=TRUE, include=FALSE}
lisas_rma <- aov_car(score ~ Error(subject / target_location * block), data = lisas)
lisas_BF <- anovaBF(score ~ target_location * block + subject, data = lisas, whichRandom = "subject")

# Block 2 ----------------------------------------------------------------------
rand <- filter(lisas, block==2 & target_location=="random")
struct <- filter(lisas, block==2 & target_location=="structured")

b2_lisas_TT <- t.test(score ~ target_location, filter(lisas, block==2), paired = TRUE)
b2_lisas_CD <- cohens_d(filter(lisas, block==2), score ~ target_location, paired = TRUE)
b2_lisas_BF <- ttestBF(rand$score, struct$score, paired = TRUE)  
```

To summarise, we observed a structured versus random condition effect in accuracy in Block 1, while in Block 2 there were no significant differences between conditions in either dependent variable. Nonetheless, across the entire experiment, and especially in Block 2, we observed nonsignificant RT effects in the direction consistent with Zhao et al. As before, to investigate behavioral performance differences in a way that is more robust to speed-accuracy tradeoffs (e.g. part of a true effect could be obscured by participants favouring accuracy over speed in Block 1, while they favoured speed over accuracy in Block 2), we calculated LISAS scores for both blocks. While we followed Zhao et al.’s original approach and report the planned t-tests for the direct replication effects in order to be able to directly compare our results, for the extended analyses we instead analysed the data using repeated measures ANOVAs to investigate potential interactions between the attentional bias effect and block. 

A repeated measures ANOVA on the LISAS scores showed strong evidence for significant main effects of block (`r apa_print(lisas_rma)$full_result[2]``r apa_print(lisas_BF)$full_result[2]`) and target location (`r apa_print(lisas_rma)$full_result[1]``r apa_print(lisas_BF)$full_result[1]`; Figure 2C). While we also found a significant interaction between block and target location (`r apa_print(lisas_rma)$full_result[3]`), the evidence for this was only anecdotal (*BF~10~* = `r as.vector(lisas_BF[4]/lisas_BF[3])`). In contrast to performance in the first block, we found no significant difference in LISAS between target locations in the second block (`r apa_print(b2_lisas_TT)$full_result`, $d$ = `r b2_lisas_CD$effsize`, *BF~10~* = `r as.vector(b2_lisas_BF)`). Overall, it therefore appears that the largest differences in LISAS between the structured and random locations occur in the first half of the experiment, with scores and performance differences between the locations stabilising by the end of Block 1 (Figure 2D).

# Participants are unaware of statistical regularities

```{r memory, message=TRUE, warning=TRUE, include=FALSE}
mem_diffs <- mem_data %>%  
  mutate(correct = ifelse(position == response, 100, 0)) %>%
  group_by(subject) %>%
  summarise(mem = mean(correct), .groups = "drop") %>% 
  right_join(diffs) %>% 
  select(subject, mem, block, rt)  
 
mem_TT <- t.test(mem ~ 1, data = filter(mem_diffs, block==1) , mu = 50) 
mem_CD <- cohens_d(mem ~ 1, data = filter(mem_diffs, block==1) , mu = 50)  
mem_BF <- ttestBF(filter(mem_diffs, block==1) $mem, mu = 50)
```

Zhao et al. found that triplets were chosen significantly more often than foils in the surprise 2AFC familiarity test (statistics from their paper: *M* = 54%, *t*(24) = 2.52, $p$ = .02, $d$ = 0.50), and concluded that this was evidence of learning of the regularities in the structured location. We did not replicate this finding - we saw no significant differences in the number of trials on which subjects chose triplets or foils as the more familiar set (`r apa_print(mem_TT)$full_result`, $d$ = `r mem_CD$effsize`, *BF~10~* = `r as.vector(mem_BF)`). Notably Zhao et al., had participants perform this test after completing the main experiment (i.e. after 450 total shape presentations), while we had them perform the test after completing the second block (i.e. after 900 shape presentations). If choosing the triplets more often on the familiarity test was indeed a reflection of learning, then we might expect that we should have an even higher proportion of trials on which triplets were chosen than reported by Zhao et al., since they would have had twice the exposure. However, our results favour the proposition that subjects did not have explicit knowledge of the regularities, as performance on the familiarity test was at chance level.

```{r memory_diffs, message=TRUE, warning=TRUE, include=FALSE}
lisas_diffs <- left_join(lisas, mem_diffs) %>%  
  spread(target_location, score) %>%
  mutate(diff = random - structured)

mem_diff_diffs <- left_join(mem_diffs, diffs) %>%  
  pivot_wider(id_cols = c(subject, mem, block), names_from = block, values_from = rt) %>% 
  mutate(rt = `2` - `1`)

mem_lisas_diffs <- left_join(lisas, mem_diffs) %>%  
  pivot_wider(id_cols = c(subject, mem, block), names_from = target_location, values_from = score) %>% 
  mutate(diff = random - structured) %>% 
  pivot_wider(id_cols = c(subject, mem), names_from = block, values_from = diff) %>% 
  mutate(diff = `2` - `1`) 

diff_b1 <- cor.test(filter(mem_diffs, block==1)$rt, filter(mem_diffs, block==1)$mem)
diff_b2 <- cor.test(filter(mem_diffs, block==2)$rt, filter(mem_diffs, block==2)$mem)

lisas_b1 <- cor.test(filter(lisas_diffs, block==1)$diff, filter(lisas_diffs, block==1)$mem)
lisas_b2 <- cor.test(filter(lisas_diffs, block==2)$diff, filter(lisas_diffs, block==2)$mem)

diff_diff <- cor.test(mem_diff_diffs$rt, mem_diff_diffs$mem)
lisas_diff <- cor.test(mem_lisas_diffs$diff, mem_lisas_diffs$mem)
```

We next wanted to see if participants who performed better on the familiarity test had a larger difference in RTs between the structured and random locations in the first block. We found no such correlation in either block (Block 1: `r apa_print(diff_b1)$full_result`; Block 2: `r apa_print(diff_b2)$full_result`). The difference of the difference between structured and random locations between the two blocks was also uncorrelated with accuracy on the familiarity test (`r apa_print(diff_diff)$full_result`), suggesting that there was no relationship between learning the regularities and changes in the attentional bias effect between blocks. Since we found attentional bias effects on both accuracy and RT, we repeated this analysis using the LISAS scores, and again found no significant correlations in either blocks (Block 1: `r apa_print(lisas_b1)$full_result`; Block 2: `r apa_print(lisas_b1)$full_result`; Figure 3A), nor for the differences between blocks (`r apa_print(lisas_diff)$full_result`; Figure 3B).

```{r proportions, message=TRUE, warning=TRUE, include=FALSE}
castor <- left_join(castor_import, distinct(select(data_raw, subject, structured_location))) %>%  
  mutate(
    accuracy_loc = as.factor(ifelse(structured_location == q_regularityloc, "Correct", "Incorrect")),
    q_regularity = as.factor(ifelse(q_regularity == 1, "Aware", "Unaware"))
    ) %>% 
  right_join(distinct(select(mem_diffs, subject, mem))) %>% 
  right_join(lisas) %>% drop_na(accuracy_loc)  

awareness_count <- castor %>%
  group_by(q_regularity) %>%
  summarise(n=n_distinct(subject))
  
loc_count <- castor %>%
  group_by(accuracy_loc) %>%
  summarise(n=n_distinct(subject))

loc_aware_count <- castor %>%
  group_by(accuracy_loc, q_regularity) %>%
  summarise(n=n_distinct(subject)) %>% 
  spread(accuracy_loc, n)

prop_accuracy <- prop.test(loc_count$n[1], sum(loc_count$n), p = 0.25) 
prop_aware_correct <- prop.test(as.matrix(loc_aware_count[2:3])) 

mem_aware_TT <- t.test(mem ~ q_regularity, castor)
mem_aware_CD <- cohens_d(castor, mem ~ q_regularity)
mem_aware_BF <- ttestBF(formula = mem ~ q_regularity, data = castor)

mem_loc_TT <- t.test(mem ~ accuracy_loc, castor)
mem_loc_CD <- cohens_d(castor, mem ~ accuracy_loc)
mem_loc_BF <- ttestBF(formula = mem ~ accuracy_loc, data = castor)
```

After completing the familiarity test, participants were asked to report whether they noticed any regularities in any of the locations. `r dc0(awareness_count$n[1]/sum(awareness_count$n)*100)`% indicated that they had noticed regularities, after which the task design was briefly explained to them. They were then asked to select which of the four locations they thought contained the structured stream (2 participants, both of whom reported not noticing regularities, did not answer this question and are excluded from the following reports). `r dc0(loc_count$n[1]/sum(loc_count$n)*100)`% of participants were able to correctly identify the location containing the structured stream, though this was not significantly different from the 25% chance level (95% CI [`r dc2(prop_accuracy$conf[1]*100)`, `r dc2(prop_accuracy$conf[2]*100)`], $p$ = `r  dc2(prop_accuracy$p.value)`, $\chi$^2^ = `r dc2(prop_accuracy$statistic[[1]])`). Participants who reported being aware of the regularities were no more likely to correctly identify their location than those who reported being unaware (aware: `r dc0(prop_aware_correct$estimate[1]*100)`%, unaware: `r dc0(prop_aware_correct$estimate[2]*100)`%; 95% CI [`r  dc2(prop_aware_correct$conf[1]*100)`, `r dc2(prop_aware_correct$conf[2]*100)`], $p$ = `r  dc2(prop_aware_correct$p.value)`, $\chi$^2^ = `r dc2(prop_aware_correct$statistic[[1]])`). There were also no differences in performance on the familiarity test between those who reported being unaware and those who reported being aware of the regularities (`r apa_print(mem_aware_TT)$full_result`, $d$ = `r mem_aware_CD$effsize`, *BF~10~* = `r as.vector(mem_aware_BF)`), nor between those who reported the location of the regularities correctly and those who did not (`r apa_print(mem_loc_TT)$full_result`, $d$ = `r mem_loc_CD$effsize`, *BF~10~* = `r as.vector(mem_loc_BF)`; Figure 3C).


```{r counterbalance, message=TRUE, warning=TRUE, include=FALSE}
# location ---------------------------------------------------------------------
lisas_counter_loc <- filter(data_raw, flag==0 & exc==0) %>%
  group_by(subject) %>%
  mutate(S_RT = sd(RT)) %>%
  group_by(subject, target_location, block, structured_location) %>%
  summarise(
    RTj = mean(RT),
    PEj = 1 - mean(correct),
    S_RT = max(S_RT),
    .groups = "drop"
  ) %>%
  group_by(subject) %>%
  mutate(S_PE = sd(PEj)) %>%
  group_by(subject, target_location, block, structured_location) %>%
  summarise(score = RTj + (S_RT / S_PE) * PEj, .groups = "drop")
 
counter_loc_rma <- aov_car(score ~ structured_location + Error(subject / target_location * block), data = lisas_counter_loc)
counter_loc_BF <- anovaBF(score ~ target_location * block * structured_location + subject, data = lisas_counter_loc, whichRandom = "subject")

# stimuli ----------------------------------------------------------------------
lisas_counter_stim <- filter(data_raw, flag==0 & exc==0) %>%
  group_by(subject) %>%
  mutate(S_RT = sd(RT)) %>%
  group_by(subject, target_location, block, structured_stream) %>%
  summarise(
    RTj = mean(RT),
    PEj = 1 - mean(correct),
    S_RT = max(S_RT),
    .groups = "drop"
  ) %>%
  group_by(subject) %>%
  mutate(S_PE = sd(PEj)) %>%
  group_by(subject, target_location, block, structured_stream) %>%
  summarise(score = RTj + (S_RT / S_PE) * PEj, .groups = "drop")
 
counter_stim_rma <- aov_car(score ~ structured_stream + Error(subject / target_location * block), data = lisas_counter_stim)
counter_stim_BF <- anovaBF(score ~ target_location * block * structured_stream + subject, data = lisas_counter_stim, whichRandom = "subject")
``` 


```{r reported_loc, message=TRUE, warning=TRUE, include=FALSE}
lisas_reported <- left_join(filter(data_raw, flag==0 & exc==0), castor) %>%
  mutate(reported_location = factor(
    ifelse(search_location == q_regularityloc, "reported structured", "reported random"),
    levels = c("reported structured", "reported random"))
    )  %>%
  group_by(subject) %>%
  mutate(S_RT = sd(RT)) %>%
  group_by(subject, reported_location, block) %>%
  summarise(
    RTj = mean(RT),
    PEj = 1 - mean(correct),
    S_RT = max(S_RT),
    .groups = "drop"
  ) %>%
  group_by(subject) %>%
  mutate(S_PE = sd(PEj)) %>%
  group_by(subject, reported_location, block) %>%
  summarise(score = RTj + (S_RT / S_PE) * PEj, .groups = "drop") %>%  
  left_join(distinct(select(castor, subject, q_regularity, accuracy_loc))) %>% 
  drop_na(accuracy_loc)

lisas_rep_loc <- aov_car(score ~ accuracy_loc + Error(subject / reported_location * block), lisas_reported) 
lisas_rep_loc_BF <- anovaBF(score ~ reported_location * block * accuracy_loc + subject, lisas_reported, whichRandom = "subject")
ls_rep_loc <- pairs(lsmeans(lisas_rep_loc, ~ reported_location | accuracy_loc | block), adjust = "bon")
 
lisas_loc <- aov_car(score ~ accuracy_loc + Error(subject / target_location * block), castor) 
ls_loc <- pairs(lsmeans(lisas_loc, ~ target_location | accuracy_loc | block), adjust = "bon")
``` 
We next considered whether the better performance in the structured location could be simply explained by a preference towards a specific location improving performance at that location, regardless of whether regularities are present. Both the location of the structured stream and the stimulus set assigned to each location were counterbalanced across participants, and we did not find any significant differences in LISAS between the locations of the structured stream (`r apa_print(counter_loc_rma)$full_result[1]`, *BF~10~* = `r as.vector(counter_loc_BF[5])`) nor between the different structured stimuli sets (`r apa_print(counter_stim_rma)$full_result[1]`, *BF~10~* = `r as.vector(counter_stim_BF[5])`). Nevertheless, if some participants had a preference for one of the random locations, their performance on the search task at that location may have benefited from a preference-driven attentional bias toward that location. However, such potential performance benefits may not be evident since the three random locations are combined into a single average. We investigated this potential scenario by examining performance on the search task based on participants' reports of where they believed (or guessed) the location of the structured stream was. That is, rather than comparisons between the *actual* structured and random locations, we instead compared performance between the location participants *thought* contained the structured stream with the others. We found significant main effects of reported location (`r apa_print(lisas_rep_loc)$full_result[2]``r apa_print(lisas_rep_loc_BF[1])$full_result`) and block  (`r apa_print(lisas_rep_loc)$full_result[4]``r apa_print(lisas_rep_loc_BF[2])$full_result`; Figure 3B). 

Crucially, when comparing performance between participants who correctly identified the structured location with participants who chose one of the random locations, we found that there was a significant interaction between the accuracy of the reported location and performance between the location participants reported contained the structured regularities and the remaining locations (`r apa_print(lisas_rep_loc)$full_result[3]``r apa_print(lisas_rep_loc_BF[10]/lisas_rep_loc_BF[6])$full_result`; Figure 3D). Neither the interaction between the reported location and block (`r apa_print(lisas_rep_loc)$full_result[6]`, *BF~10~* = `r as.vector(lisas_rep_loc_BF[4]/lisas_rep_loc_BF[3])`), nor the three-way interaction between accuracy, reported location and block were significant (`r apa_print(lisas_rep_loc)$full_result[7]`, *BF~10~* = `r as.vector(lisas_rep_loc_BF[18]/lisas_rep_loc_BF[17])`). 

Follow-up Bonferroni-adjusted pairwise tests showed that for participants who accurately identified the structured location, LISAS scores were significantly lower in the structured location than in the random locations in both the first (`r apa_print(ls_rep_loc)$full_result[1]`) and second (`r apa_print(ls_rep_loc)$full_result[3]`) blocks. For participants who indicated that one of the random locations contained the regularities, performance in that chosen (random) location was only significantly better than their performance in the other locations in the second block (`r apa_print(ls_rep_loc)$full_result[4]`), and not in the first (`r apa_print(ls_rep_loc)$full_result[2]`). When repeating this analysis but instead looking at performance differences between the structured and random locations (Figure 3D), we found that in both those who correctly reported the structured location (`r apa_print(ls_loc)$full_result[1]`) and those who chose one of the random locations (`r apa_print(ls_loc)$full_result[2]`) had significantly lower LISAS in the structured location compared to the random locations in the first block. In the second block, only those who correctly reported the structured location showed a significant difference (`r apa_print(ls_loc)$full_result[3]`), while those who reported that the regularities occurred in a random location showed no significant performance differences between the structured and random locations (`r apa_print(ls_loc)$full_result[4]`). Together, these results indicate that it is indeed the regularities and not a location-based preference that biases attention towards the structured location in the first half of the experiment. In the second half, the attentional bias effect diminishes but persists for some participants, while others shift their attention towards one of the random locations, perhaps as they search to uncover new or more complex structure in the environment.


# Discussion
We conducted a direct replication of Zhao et al.’s [@Zhao2013] Experiment 1, more than tripling the sample size from the original N = 25 to ensure sufficient statistical power to detect subtle effects. We expected that the results would provide an insight into the strength of evidence behind the claim that attention is automatically biased towards task-irrelevant regularities in the environment. While we could not directly replicate Zhao et al.’s finding that response times are faster to visual search targets presented in a location containing statistical regularities, we did conceptually replicate their effect. Our results also show that performance on the search task was better in the location containing regularities, though we found that this was primarily reflected in accuracy rates rather than response speed. A potential reason for why we found the facilitation effect in accuracy rather than response speed may be as simple as subtle differences in task instructions, which have been previously shown to modulate speed-accuracy trade-offs [@vanVeen2008]. The main result of our replication therefore supports their conclusion that attention is spontaneously biased to locations that contain statistical regularities.
 
To extend Zhao et al.’s initial findings, here we investigated whether the magnitude of the attention bias effect decreases with exposure to the regularities. In the second half of the experiment, we did not see a difference in performance between the two target locations, neither in accuracy nor RTs. When accounting for potential speed-accuracy trade-offs by integrating both RTs and accuracy into a single weighted score, our results show that overall performance was better in the structured relative to the random locations at the beginning of the experiment, with performance differences between conditions decreasing over time. This result is consistent with studies showing that attentional biases towards regularities occur at an early stage of exposure, gradually reducing over time [@Alamia2016; @Yu2015], suggesting that the attentional bias towards regularities diminishes with continued exposure. It is possible that attention fluctuates towards learning opportunities in the environment, leading to a diminishing attentional bias towards regularities as they are learned and expectations in that location are strengthened. It remains unclear whether continued additional exposure would entirely eliminate all performance differences between conditions. Perhaps once the regularities in the structured location are learned we may even see a reversal, with attention being biased away in the search for something new to learn. In any case, this finding demonstrates the critical need for research into the relationship between attention and statistical learning to consider the temporal dynamics of their interplay. Lack of such consideration may also help to explain some of the mixed results reported in the literature. The potential nuances of this relationship would be best investigated by future studies that manipulate the amount of learning opportunities and directly probe individual learning rates with online or periodic measures.

Zhao et al. additionally reported evidence of explicit knowledge of the regularities in the structured location: in a surprise 2AFC familiarity test conducted after the main task, triplets of shapes from the structured location were selected more frequently than foils. We did not replicate this finding, and instead found substantial evidence for the hypothesis that learning was implicit, as participants could not discriminate between triplets and foils above chance levels. This failure to replicate is somewhat surprising considering that participants had twice as much exposure to the regularities in our replication than in Zhao et al.’s original study. Given that we also found evidence in support of an attentional bias towards the structured location, we would expect accuracy on the familiarity test to be higher than that reported by Zhao et al. due to the additional exposure. 

One possible explanation for this discrepancy is perhaps that as the attentional bias towards regularities diminishes with exposure, the regularities participants might have learned at the start of the experiment are simply forgotten once attention starts to shift away; perhaps if we had conducted the familiarity test after the first block (matching Zhao et al.) we would have also found evidence of explicit learning. However, our results don’t provide support for this explanation, as we found no relationship between accuracy on the familiarity test and the strength of the attentional bias effect on performance in either block, indicating that participants who exhibited larger biases towards regularities in the second half of the experiment were no more likely to show explicit knowledge of the regularities than those participants who no longer showed a bias towards regularities in the second block. The lack of learning we found here is however consistent with previous studies that have demonstrated that only attended regularities are learned [@Alamia2016; @Baker2004; @Toro2005; @Damsma2020], and suggests that despite attention being biased towards regularities, selective attention may be necessary for statistical learning to occur. 

One limitation of both the original study and our replication is that strong conclusions as to the automaticity and implicitness of the attentional bias cannot be made. Since the sequences were simply viewed passively in between search trials, intentional processing of the structure cannot be ruled out. Although participants did not explicitly learn the regularities nor their location, as evidenced by the chance performance on both the familiarity test and in reporting the location of the structured location, it remains possible that observers may have nonetheless directly attended to the streams in an attempt to figure out the underlying structure. Such direct attention may have been sufficient to produce the attentional bias effect, despite not being enough for the regularities to be learned to an extent that would be captured by our measures.  

Zhao et al. argued that their results support the attentional bias being an implicit and automatic process, rather than participants simply becoming explicitly aware of the regularities and intentionally orienting to them. In support of this argument, they highlight that although performance was above chance on the surprise familiarity test, only 7 participants (30%) reported noticing the regularities, of whom only 1 (14%) also correctly identified the structured location. While we found chance-level performance on the familiarity test, we also found no evidence for explicit awareness of the regularities, as even participants who reported being aware of regularities did not identify the structured location above chance levels. However, we found that the attentional bias towards regularities was strongest in the first half of the experiment, and only persisted for those participants who correctly identified the structured location. It is therefore possible that at the beginning of the experiment, attention was automatically biased towards the regularities, and subsequently some participants continued to direct their attention to that location due to some implicit awareness of the regularities contained within that location. Alternatively, it is possible that rather than an awareness of the regularities, participants had some awareness of the attentional bias towards one of the locations, and reported the regularities in whichever location they were attending to more towards the end of the experiment. Our results therefore do not rule out the possibility that the attentional bias effect was at least partially driven by awareness and top-down attentional selection of the location containing regularities. However, they do suggest that awareness of the regularities is not necessary for attention to be biased towards them. Future studies will need to investigate this question more directly before any strong conclusions regarding the level of awareness and intentionality can be drawn.

## Conclusion
In summary, while our high-powered replication attempt did not directly reproduce Zhao et al.’s central finding that response times are faster to visual search targets presented in a location containing task-irrelevant regularities, our results offer a clear conceptual replication. We found improved performance (mostly in accuracy) in the locations containing regularities, supporting their original conclusion that attention is biased towards task-irrelevant statistical regularities in the environment. We extend upon Zhao et al.’s original work by doubling the duration of the experiment to investigate whether the attentional bias effect diminishes with continued exposure. Our results indicate that this is indeed the case, as we found that while performance was initially improved in the structured location, these performance benefits dissipated over time and almost disappeared by the end of the experiment. Despite doubling the amount of exposure to the regularities, we failed to replicate Zhao et al.’s finding that these regularities were explicitly learned. Together, our results suggest that while initially attention may be spontaneously bias towards regularities in the environment, this attentional bias may be transient. It remains unclear whether this is due to rapid learning of the contingencies early on, after which the environment continues to be explored in the search for new learning opportunities, or whether it is due to attention drifting if the contingencies are too difficult to learn and are otherwise not task-relevant. Future research could provide new insights by investigating how the amount of learning and learning opportunities in the environment may interactively influence the attentional bias effect over time.

<!-- In all cases, observers had seen the individual shapes in the triplet versus foil sequences an equal number of times, and so the sequences were only distinguished by the familiarity of their higher-order temporal structure. -->

### Acknowledgements
We thank Victoria Suess for her assistance with data acquisition. 

### Ethical Statement
The experiments were approved by the local ethics committee (CMO Arnhem-Nijmegen, The Netherlands) under the general ethics approval (CMO 2014/288, v. 2.1) and were conducted in compliance with these guidelines. All participants provided written informed consent.

### Funding Statement
This study was supported by is supported by the Netherlands Organisation for Scientific Research (Vidi grant 406.16.535) and the European Research Council (ERC-StG 678286 and ERC-CoG 101000942) awarded to F.P.d.L, The Netherlands Organisation for Scientific Research (NWO Veni grant 016.Veni.198.065) awarded to E.S., and the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement (H2020-MSCA-IF-2018 No. 846702), awarded to A.V.

### Data accessibility
All data and code used for stimulus presentation and analysis are released publicly under an open license at the Donders Repository: https://data.donders.ru.nl/collections/di/dccn/DSC_3018028.05_239.

### Competing Interests
We declare we have no competing interests.

### Author's Contributions
All authors designed research. A.V. performed research, analysed data, drafted the manuscript and made figures. All authors edited and approved the final version of the manuscript.

# Figures

```{r plotting, include=FALSE}
fillcols <- c("#EE6677","#4477AA")

# Common theme -----------------------------------------------------------------
mytheme <- theme_cowplot() +
  theme(
    axis.text = element_text(size = 10),
    axis.title = element_text(size = 12),
    axis.ticks.length = unit(0.1, "cm"),
    legend.key.size = unit(.4, "cm"),
    legend.text = element_text(size = 10),
    legend.position = "top",
    legend.direction = "horizontal",
    legend.justification = "center"
  ) +
  background_grid(major = "y", minor = "y", size.major = .1, size.minor = .1)
 
# Common geoms -----------------------------------------------------------------
violin <- geom_violin(
  trim = FALSE, scale = "count",
  alpha = 1, linetype = 0,
  position = position_dodge(width = .9), width = 1
)

box <- geom_boxplot(
  aes(group=interaction(block,target_location)),
  position = position_dodge(width = .9), width = .08,
  fill = "white",
  alpha = 1,
  notch = TRUE, notchwidth = .2,
  outlier.size = .25, weight = .5
)
```

```{r figure1, echo=FALSE, fig.cap="Task design. *A. Stimuli.* Streams of shapes were shown at the top, bottom, left and right of a central fixation cross, with each stream generated from a different set of nine shapes. One of the locations contained temporal regularities (structured location, shown in red); the shapes were grouped into three triplets, with the shapes in each triplet always presented in the same order, and the order of the triplets shuffled. In the remaining three locations the shapes were presented in a random order (random location, shown in blue). *B. Timeline.* The shape streams were periodically interrupted by a visual search array containing one target 'T' shape and three 'L'-shaped distractors. The task was to indicate whether the target was pointing leftwards or rightwards (the target in the example illustrated here is pointing leftwards). After a response was made, the shape streams would resume until the next visual search trial."}

# Fig 1A: Timeline ------------------------------------------------------------- 
p1 <- readPNG("figures/timeline.png")
g1 <- grid::rasterGrob(p1)
timeline <- ggplot() +
  annotation_custom(g1) +
  labs(tag = "B") +
  mytheme + theme(axis.line=element_blank())

# Fig 1B: Stimulus ------------------------------------------------------------- 
p2 <- readPNG("figures/stimulus.png")
g2 <- grid::rasterGrob(p2, interpolate = TRUE,
                       x = unit(.5, "npc"), y = unit(.5, "npc"))
stimulus <- ggplot() +
  annotation_custom(g2) +
  labs(tag = "A") +
  mytheme + theme(axis.line=element_blank())

# Combine & save --------------------------------------------------------------- 
png("figures/figure1.png", width = 15, height = 8, units = "cm", res = 300)
stimulus + timeline + plot_layout(widths = c(.4,.6))
dev.off()
include_graphics("figures/figure1.png")
```

```{r figure2, echo=FALSE, fig.cap="*A. Response times on the visual search task.* Mean response times on the visual search task are shown separately for Block 1 (left) and 2 (right), split by whether the target was presented in the structured (red) or random (blue) condition. The boxplot indicates the median response time, with the notch highlighting the 95% confidence bounds and the boundary indicating the top and bottom quartiles. Whiskers indicate the 10th and 90th percentiles and the points above and below the whiskers show outliers who fall outside this range. *B. Accuracy on the visual search task.* Mean percent correct on the visual search task is shown for Block 1 (left) and 2 (right), split by whether the target was presented in the structured (red) or random (blue) condition. The boxplot indicates the median percent correct; other plotting details are as indicated directly above. *C. Linear integrated speed accuracy (LISAS) on the visual search task.*Mean LISAS on the visual search task are shown for Block 1 (left) and 2 (right), split by whether the target was presented in the structured (red) or random (blue) condition. The boxplot indicates the median LISAS; other plotting details are as indicated directly above. *D. Mean LISAS as a function of trial.* Trials were binned into 9 bins (20 trials each), and the mean LISAS for the visual search task was calculated for each bin. Points show mean LISAS for each bin for both the structured (red) and random (blue) conditions, with error bars indicating +/-SEM."}

# Fig 2A: RT ------------------------------------------------------------------- 
g.rt <-
  ggplot(data_means, aes(y = rt, x = block, fill = target_location)) +
  mytheme + violin + box + 
  scale_color_manual(values = fillcols, guide = NULL) +
  scale_fill_manual(
    values = fillcols,
    guide = guide_legend(override.aes = list(linetype = 0, alpha = 1, size = 1))
  ) +
  scale_y_continuous(breaks = seq(.5, 3, .5), labels = sprintf("%.1f", seq(.5, 3, .5))) +
  scale_x_discrete(labels = c("1" = "Block 1", "2" = "Block 2")) +
  labs(x = "", y = "Response Time (s)", fill = "", tag = "A") + 
  geom_signif(
    y_position = 3.16, xmin = .775, xmax = 1.225, annotation = paste(""),
    tip_length = c(.06, .06), size = .25
  ) +
  geom_signif(
    y_position = 3.16, xmin = 1.775, xmax = 2.225, annotation = paste(""),
    tip_length = c(.06, .06), size = .25
  ) +
  geom_signif(
    y_position = 3.33, xmin = 1, xmax = 2, annotation = paste("*"),
    tip_length = c(.075, .075), size = .25, vjust = .5
  ) 

# Fig 2B: Accuracy -------------------------------------------------------------  
g.acc <- 
  ggplot(data_means, aes(y = accuracy, x = block, fill = target_location)) +
  mytheme + violin + box + 
  scale_color_manual(values = fillcols, guide = NULL) +
  scale_fill_manual(values = fillcols, guide = NULL) +
  scale_y_continuous(breaks = seq(40, 100, 10), labels = sprintf("%.0f", seq(40, 100, 10))) +
  scale_x_discrete(labels = c("1" = "Block 1", "2" = "Block 2")) + 
  labs(x = "", y = "Percent correct", fill = "", tag = "B") + 
  coord_cartesian(ylim = c(32,125), expand = FALSE) +
  geom_signif(
    y_position = 116, xmin = .775, xmax = 1.225, annotation = paste("***"),
    tip_length = c(.06, .06), size = .25, vjust = 1.5
  ) +
  geom_signif(
    y_position = 116, xmin = 1.775, xmax = 2.225, annotation = paste(""),
    tip_length = c(.06, .06), size = .25
  ) +
  geom_signif(
    y_position = 120.4, xmin = 1, xmax = 2, annotation = paste("*"),
    tip_length = c(.075, .075), size = .25, vjust = .5
  ) 
 
# Fig 2C: LISAS -------------------------------------------------------------------
g.lisas <-
  ggplot(lisas, aes(y = score, x = block, fill = target_location)) + mytheme +
  geom_violin(
    trim = FALSE, scale = "count",
    alpha = 1, linetype = 0,
    position = position_dodge(width = .9), width = .8
  ) + 
  box +  
  scale_color_manual(values = fillcols, guide = NULL) +
  scale_fill_manual(values = fillcols, guide = NULL) +
  scale_y_continuous(breaks = seq(0, 7, 1), labels = sprintf("%.1f", seq(0, 7, 1))) +
  scale_x_discrete(labels = c("1" = "Block 1", "2" = "Block 2")) +
  labs(x = "", y = "LISAS", fill = "", tag = "C") +
  geom_signif(
    y_position = 6.94, xmin = .775, xmax = 1.225, annotation = paste("***"),
    tip_length = c(.06, .06), size = .25, vjust = 1.5
    ) +
  geom_signif(
    y_position = 6.94, xmin = 1.775, xmax = 2.225, annotation = paste(""),
    tip_length = c(.06, .06), size = .25
    ) +
  geom_signif(
    y_position = 7.3, xmin = 1, xmax = 2, annotation = paste("*"),
    tip_length = c(.07, .07), size = .25, vjust = .5
    )

# Fig 2D: Binned LISAS ---------------------------------------------------------
lisas_binned <- filter(data_raw, flag==0 & exc==0) %>%
  mutate(bin = cut(trial, breaks = seq(0, 160, 20), labels = FALSE)) %>% 
  group_by(subject) %>%
  mutate(S_RT = sd(RT)) %>%
  group_by(subject, target_location, bin) %>%
  summarise(
    RTj = mean(RT),
    PEj = 1 - mean(correct),
    S_RT = max(S_RT),
    .groups = "drop"
  ) %>%
  group_by(subject) %>%
  mutate(S_PE = sd(PEj)) %>%
  group_by(subject, target_location, bin) %>%
  summarise(score = RTj + (S_RT / S_PE) * PEj, .groups = "drop") %>% 
  group_by(target_location, bin) %>%
  summarise(mscore = mean(score), sem = sd(score) / sqrt(n()), .groups = "drop")

g.lisas_bin <- ggplot(lisas_binned, aes()) + mytheme + 
  geom_vline(xintercept = 4.5, linetype = "dashed", colour = "darkgrey") + 
  geom_line(data = lisas_binned, aes(bin, mscore, color = target_location))+
  geom_point(data = lisas_binned, aes(bin, mscore, color = target_location), size = 1) +
  geom_errorbar(data = lisas_binned,
                aes(x = bin, y = mscore, color = target_location,
                    ymin = mscore - sem / 2, ymax = mscore + sem / 2), width = .2) +
  scale_colour_manual(values = fillcols, guide = NULL) +
  scale_fill_manual(values = fillcols, guide = NULL) +
  labs(y = "LISAS", colour = "", tag = "D") +
  scale_x_continuous(
    name = "Trial number",
    breaks = seq(1, 8, 1),
    labels = c(sprintf("%.0f", seq(20, 160, 20)))
    )  

# Combine & save ---------------------------------------------------------------
png("figures/figure2.png", width = 18, height = 16, units = "cm", res = 300)
g.rt + g.acc + g.lisas + g.lisas_bin +
  plot_layout(
  design = "
  AB
  CD
  ",
  guides = "collect") &
  theme(legend.position = "bottom", legend.direction = "horizontal")
dev.off()
include_graphics("figures/figure2.png")
```

```{r figure3, echo=FALSE, fig.cap="*A. Correlation between accuracy on the Familiarity test and the difference in LISAS between conditions in Block 1.* The y-axis shows the difference in LISAS between conditions in Block 1 for each participant (grey crosses), calculated by subtracting each participant’s mean LISAS in the structured location from the mean LISAS in the random locations. Positive values indicate better performance in the structured location than in the random. The difference score is plotted against each participant’s percent correct score on the Familiarity test on the x-axis. The linear regression line is shown in black. *B. Correlation between accuracy on the Familiarity test and the change between blocks in the difference in LISAS between conditions.* The y-axis shows the change between blocks in the LISAS difference between conditions for each participant (grey crosses), calculated by subtracting the difference scores plotted in panel C (Block 1) from each participant’s difference scores in Block 2. Positive values indicate that performance gains in the structured location increased between blocks. The difference score is plotted against each participant’s percent correct score on the Familiarity test on the x-axis. The linear regression line is shown in black. *C. Accuracy on the Familiarity test split by accuracy of reported structured location.* Mean percent correct is shown separately for subjects who reported the location of the structured stream correctly (n = 25) and for those who didn’t (n = 50). The notches highlight the 95% confidence bounds for the median accuracy, the boundaries indicate the top and bottom quartiles, and the whiskers indicate the 10th and 90th percentiles. The dashed horizontal line shows the chance level of performance (50% correct). D. Performance on the visual search task in the reported locations, split by accuracy of the report. Mean LISAS scores for the location reported as structured (red) and the locations reported as random (blue) are plotted separately for participants who correctly identified the structured location vs those who did not. The notch indicates the median LISAS, with the notch highlighting the 95% confidence bounds and the boundary indicating the top and bottom quartiles. Whiskers indicate the 10th and 90th percentiles and the points above the whiskers show outliers who fall outside this range."}

# Fig 3C: LISAS (rand-struct) --------------------------------------------------  
g.scat <- lisas_diffs %>% 
  ggscatter(x = "mem", y = "diff", size = 1, shape = 4, color = "grey50",
            add = "reg.line", add.params = list(size=.8, color = "black")) +
  stat_cor( label.x.npc = .65, label.y.npc = 1, size = 2.5) +
  labs(x = "Percent correct", y = bquote(Delta * " LISAS"), tag = "A") +
  mytheme + background_grid(major = "none", minor = "none")  

# Fig 3D: B2(rand-struct) - B1(rand-struct) ------------------------------------ 
g.scat2 <- cbind(lisas_diffs, mem_lisas_diffs$diff) %>%  
  rename(diffdiff = "mem_lisas_diffs$diff") %>% 
  ggscatter(x = "mem", y = "diffdiff", size = 1, shape = 4, color = "grey50", 
            add = "reg.line", add.params = list(size=.8, color = "black")) +
  stat_cor( label.x.npc = .65, label.y.npc = 1, size = 2.5) +
  labs(x = "Percent correct", 
       y = bquote("B2"[Delta * "LISAS "] * - " B1"[Delta * "LISAS"]), tag = "B") +
  mytheme + background_grid(major = "none", minor = "none") + 
  theme(axis.title.y = element_text(size = 10),)

# Fig 3A: Location accuracy ---------------------------------------------------- 
g.loc_mem <- filter(castor, block==1) %>%   
  ggplot(aes(accuracy_loc, mem)) + mytheme +
  geom_hline(yintercept = 50, linetype = "dashed") +
  geom_violin(fill = "grey", trim = FALSE, scale = "count", alpha = .7, linetype = 0) +
  geom_boxplot(fill = "white", width = .2, notch = TRUE, notchwidth = .2, weight = .5) +
  labs(x = "Reported structured location", y = "Percent correct", tag = "C")  

# Fig 3D: LISAS x location accuracy --------------------------------------------  
g.sub_lisas <- lisas_reported %>% 
  group_by(subject, reported_location, accuracy_loc) %>%
  summarise(mscore = mean(score), .groups = "drop") %>%
  ggplot(aes(x = accuracy_loc, y = mscore, fill = reported_location)) + 
  mytheme + violin +
  geom_boxplot(
    aes(group = interaction(accuracy_loc, reported_location)),
    position = position_dodge(width = .9), width = .08,
    notch = TRUE, notchwidth = .2, outlier.size = .25, weight = .5,
    show.legend = FALSE, fill = "white"
  ) +
  scale_fill_manual(
    values = fillcols,
    guide = guide_legend(override.aes = list(linetype = 0, alpha = 1, size = 1))
  ) +
  labs(x = "Reported structured location", y = "LISAS", tag = "D", fill = "") +
  geom_signif(
    y_position = 6, xmin = .775, xmax = 1.225, annotation = paste("***"),
    tip_length = c(.06, .06), size = .25, vjust = 1.5
  ) +  
  geom_signif(
    y_position = 6.42, xmin = 1, xmax = 2, annotation = paste("*"),
    tip_length = c(.1,.1), size = .25, vjust = .5
  )

# Combine & save ---------------------------------------------------------------
png("figures/figure3.png", width = 16, height = 15, units = "cm", res = 300)
g.scat + g.scat2 + (ggplot() + theme_void()) + g.loc_mem + g.sub_lisas +
  plot_layout(design = c(
    area(1,1,8,4),
    area(1,5,8,9),
    area(9,1,9,9),
    area(10,1,22,3),
    area(10,4,22,9) 
  ))
dev.off()
include_graphics("figures/figure3.png")

```
 
 

# References
\printbibliography 
