Our Experiment: Each eyetrackingR vignette uses theeyetrackingRpackage to analyze real data from a simple 2-alternative forced choice (2AFC) word recognition task administered to 19- and 24-month-olds. On each trial, infants were shown a picture of an animate object (e.g., a horse) and an inanimate object (e.g., a spoon). After inspecting the images, they disappeared and they heard a label referring to one of them (e.g., “The horse is nearby!”). Finally, the objects re-appeared on the screen and they were prompted to look at the target (e.g., “Look at the horse!”).

In this vignette, we want to examine our data *without* throwing out time as a predictor. That is, we want to ask how are differences (if any) between the Target conditions emerged over time in each trial. To do so, we will perform a growth curve analysis.

Before performing this analysis, we’ll need to prepare and clean our dataset. Here we will to do this quickly and with few notes but, for more information, see the vignette on preparing your data.

```
library("Matrix")
library("lme4")
library("ggplot2")
library("eyetrackingR")
data("word_recognition")
data <- make_eyetrackingr_data(word_recognition,
participant_column = "ParticipantName",
trial_column = "Trial",
time_column = "TimeFromTrialOnset",
trackloss_column = "TrackLoss",
aoi_columns = c('Animate','Inanimate'),
treat_non_aoi_looks_as_missing = TRUE
)
# subset to response window post word-onset
response_window <- subset_by_window(data,
window_start_time = 15500,
window_end_time = 21000,
rezero = FALSE)
```

`## Avg. window length in new data will be 5500`

```
# analyze amount of trackloss by subjects and trials
(trackloss <- trackloss_analysis(data = response_window))
# remove trials with > 25% of trackloss
response_window_clean <- clean_by_trackloss(data = response_window, trial_prop_thresh = .25)
```

`## Performing Trackloss Analysis...`

`## Will exclude trials whose trackloss proportion is greater than : 0.25`

`## ...removed 33 trials.`

```
# create Target condition column
response_window_clean$Target <- as.factor( ifelse(test = grepl('(Spoon|Bottle)', response_window_clean$Trial),
yes = 'Inanimate',
no = 'Animate') )
```

In this analysis, we are interested in the change of our data over time. Therefore, it’s very important that we visualize our data to see whether our statistical estimates make any sense at all.

To plot the time-course of our data, we first need to convert it to eyetrackingR’s `time_sequence_data`

format, which summarizes the data into time-bins and calculates proportion-looking for each (here we use 100ms time bins). From there, we simply use the `plot`

method to plot the time-course.

```
# aggregate across trials within subjects in time analysis
response_time <- make_time_sequence_data(response_window_clean, time_bin_size = 100,
predictor_columns = c("Target"),
aois = "Animate"
)
# visualize time results
plot(response_time, predictor_column = "Target") +
theme_light() +
coord_cartesian(ylim = c(0,1))
```

From this plot, we should expect to see some big differences between Target conditions over time. However, an ideal analysis will also be consistent with the observation that these differences emerged after ~500ms – i.e, they were not present at window onset.

Growth curve analysis (GCA) lets us model the timecourse of attention by fitting curves to proportion-looking over the course of the trial, and statistically assessing the bends in these curves. Our implementation of a growth curve analysis modelled after Mirman et al. (2008).

eyetrackingR sets us up nicely for GCA. Above, we used `make_time_sequence_data`

to generate the dataframe. This dataframe includes everything we need for GCA.

First, it includes the same dependent variable columns as `make_time_window_data`

, giving us an option of analyzing raw proportions or its transformations (for more detail, see the documentation or the window analysis vignette ).

`Prop`

– the mean of raw proportion scores`LogitAdjusted`

– the logit transformation,`log( Prop / (1-Prop) )`

, adjusted to avoid -/+ infinity`Elog`

– the empirical logit transformation`log( Prop+e / (1-Prop+e) )`

`ArcSin`

– the arcsine-root transformation

Second, `time_sequence_data`

has a series columns corresponding to ‘orthogonal polynomial timecodes.’ You can think of these as the linear, quadratic, cubic, etc. component of our `Time`

predictor (so `Time`

, `Time^2`

, `Time^3`

, etc.). However, unlike simply taking the power, these transformations are *uncorrelated* with each other, and therefore more appropriate for multiple regression.

If this seems confusing, it might help to simply visualize each of these vectors:

```
# generate dataframe summarizing values of each vector
timecodes <- unique(response_time[, c('ot1','ot2','ot3','ot4','ot5','ot6','ot7')])
timecodes$num <- 1:nrow(timecodes)
ggplot(timecodes, aes(x=num, y=ot1)) +
geom_line() +
geom_line(aes(y=ot2), color='red') + # quadratic
geom_line(aes(y=ot3), color='blue') + # cubic
geom_line(aes(y=ot4), color='green') + # quartic
geom_line(aes(y=ot5), color='purple') + # quintic
geom_line(aes(y=ot6), color='yellow') + # sextic
geom_line(aes(y=ot7), color='pink') + # septic
scale_x_continuous(name="") +
scale_y_continuous(name="")
```

`round(cor(timecodes[, c(1:7)]),5)`

```
## ot1 ot2 ot3 ot4 ot5 ot6 ot7
## ot1 1 0 0 0 0 0 0
## ot2 0 1 0 0 0 0 0
## ot3 0 0 1 0 0 0 0
## ot4 0 0 0 1 0 0 0
## ot5 0 0 0 0 1 0 0
## ot6 0 0 0 0 0 1 0
## ot7 0 0 0 0 0 0 1
```

You can see the linear time code moving from the bottom left to top right corner of that plot (in black). Each of the other coloured lines corresponds to a different growth function.

The idea behind GCA is that we can simulanteously regress differences between our conditions (here, Target) on each of these to see which (in combination, or independently) best captures the pattern of growth in our data. Importantly, because they are uncorrelated, they are sure to capture *distinct* variance in our data.

Let’s fit our first GCA model, **including only the linear time code for now** (i.e., we temporarily want to ignore non-linear change over time):

```
# sum-code and center our predictor:
response_time$TargetC <- ifelse(response_time$Target == 'Animate', .5, -.5)
response_time$TargetC <- as.numeric(scale(response_time$TargetC, center=TRUE, scale=FALSE))
# Construct model
model_time_sequence <- lmer(Elog ~ TargetC*(ot1) + (1 + ot1 | Trial) + (1 + ot1 | ParticipantName),
data = response_time, REML = FALSE)
# cleanly show important parts of model (see `summary()` for more)
broom::tidy(model_time_sequence, effects = "fixed")
```

```
## term estimate std.error statistic
## 1 (Intercept) 0.7402578 0.1422129 5.205280
## 2 TargetC 1.7601735 0.1388628 12.675630
## 3 ot1 -2.9962321 0.6338293 -4.727191
## 4 TargetC:ot1 -0.9950565 0.7320837 -1.359211
```

`drop1(model_time_sequence, ~., test="Chi")`

```
## Single term deletions
##
## Model:
## Elog ~ TargetC * (ot1) + (1 + ot1 | Trial) + (1 + ot1 | ParticipantName)
## Df AIC LRT Pr(Chi)
## <none> 28106
## TargetC 1 28123 19.4738 1.02e-05 ***
## ot1 1 28118 14.0063 0.0001822 ***
## TargetC:ot1 1 28105 1.6244 0.2024827
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

eyetrackingR includes a `plot`

method for `time_sequence_data`

that can easily overlay the predictions of this model on the raw data:

```
plot(response_time, predictor_column = "Target", dv = "Elog", model = model_time_sequence) +
theme_light()
```

This model reveals two main effects of TargetC and ot1, respectively. These may seem OK at first, but they are actually not a great match to our data (as revealed by the plot of their predictions).

- The main effect of TargetC says that, on average throughout each trial, participants looked more to the Animate target when it was named than when the Inanimate was named. This is accurate.
- The main effect of ot1 says that, as the trial progressed, participants tended to look away from the Animate over the course of each trial. This is also accurate.

But what we are missing is a result that reflects that these differences between Target *emerged* after the start of the window.

To remedy this issue, we can allow for non-linear change overtime by including a few higher-order orthogonal polynomials:

```
model_time_sequence <- lmer(Elog ~ TargetC*(ot1 + ot2 + ot3 + ot4) + (1 | Trial) + (1 | ParticipantName),
data = response_time, REML = FALSE)
# commented out because this many random slopes takes too long to fit for a vignette
# model_time_sequence <- lmer(Elog ~ TargetC*(ot1 + ot2 + ot3 + ot4) +
# (1 + ot1 + ot2 + ot3 + ot4 | Trial) + (1 + ot1 + ot2 + ot3 + ot4 | ParticipantName),
# data = response_time, REML = FALSE)
# cleanly show important parts of model (see `summary()` for more)
broom::tidy(model_time_sequence, effects = "fixed")
```

```
## term estimate std.error statistic
## 1 (Intercept) 0.7336737 0.1429259 5.133245
## 2 TargetC 1.7553628 0.1383939 12.683814
## 3 ot1 -2.8241968 0.1930699 -14.627847
## 4 ot2 0.3190222 0.1933733 1.649774
## 5 ot3 -0.5844908 0.1928014 -3.031569
## 6 ot4 -0.3245295 0.1921740 -1.688727
## 7 TargetC:ot1 -0.7895461 0.4012825 -1.967557
## 8 TargetC:ot2 -1.7446826 0.4020504 -4.339462
## 9 TargetC:ot3 4.2166751 0.4009058 10.517871
## 10 TargetC:ot4 -1.8162548 0.3998975 -4.541800
```

`drop1(model_time_sequence, ~., test="Chi")`

```
## Single term deletions
##
## Model:
## Elog ~ TargetC * (ot1 + ot2 + ot3 + ot4) + (1 | Trial) + (1 |
## ParticipantName)
## Df AIC LRT Pr(Chi)
## <none> 28096
## TargetC 1 28114 19.327 1.101e-05 ***
## ot1 1 28305 210.511 < 2.2e-16 ***
## ot2 1 28097 2.721 0.099027 .
## ot3 1 28104 9.184 0.002442 **
## ot4 1 28097 2.851 0.091310 .
## TargetC:ot1 1 28098 3.870 0.049159 *
## TargetC:ot2 1 28113 18.804 1.449e-05 ***
## TargetC:ot3 1 28204 109.693 < 2.2e-16 ***
## TargetC:ot4 1 28115 20.595 5.674e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

```
plot(response_time, predictor_column = "Target", dv = "Elog", model = model_time_sequence) +
theme_light()
```

This model better captures our data because it shows critical interactions between Target and higher-order polynomials. These allow the model to capture the steep drop in the Inanimate condition as well as the bend that follows it.

GCA can tell you that your conditions/predictors had some predictive value (i.e., changed differently) over time. They can also tell you the form of that change by seeing the polynomials with which they reliably interacted. However, they can’t tell you at what time these predictors had an effect.

In the case of ascertaining reaction times, one way to answer this question is by using onset-contingent analyses.

For more general approaches to estimating time windows of divergence (estimating both the onset and offset of divergences), take a look at the divergence vignette.

Mirman, D., Dixon, J., & Magnuson, J. S. (2008). Statistical and computational models of the visual world paradigm: Growth curves and individual differences. Journal of Memory and Language, 59, 474–494.