Because the gaze samples tend to be noisy, different methods use
various approaches for computing velocity from noisy samples. Methods by
Engbert and Kliegl (2003) and Otero-Millan et al. (2014) used the same
approach based on averaging over multiple samples to compute velocity,
whereas Nyström and Holmqvist (2010)
compute a simple derivative and then filter it. By default, package uses
the former approach (`velocity_function = diff_ek`

) but you
can also use the latter (`velocity_function = diff_nh`

) or
implement a custom method (see below). Acceleration is computed the same
way but from velocity samples.

This method computes velocity over a predefined time window
controlled via `ek_velocity_time_window`

that defaults to 20
ms. The width of the window is converted to samples based on the
`sampling_rate`

parameter and is ensured to be an
*odd* number equal to or larger than 3. I.e., for a window that
is three samples long, the velocity is computed based on one sample
before and one sample after the current, for a five-samples window - two
samples before and after, etc. The code always starts with the requested
window size but all iteratively reduces it by two samples to accommodate
for trial limits and missing values. An `NA`

velocity is
returned only if even the smallest window (three samples wide) cannot be
used.

In general, horizontal and vertical velocity components are computed
as \[v_x[i] =
\frac{\sum_{j=1}^{(N-1)/2}x[i+j] -
x[i-j]}{\sum_{j=1}^{(N-1)/2}2j\cdot\Delta t}\] where \(i\) is the index of a sample, \(\Delta t = \frac{1}{sample~rate}\) is a
duration of a single sampling frame, and \(N\) is an *odd* integer width of the
moving average used to compute the velocity. In Engbert and Kliegl (2003), \(N=5\) and \(\Delta t = 4\) (250 Hz sampling rate)
translate into a 20 ms moving average window (default value used in the
method implementation). Below is a derivation that shows the equivalence
of the formula above to formula 1 in Engbert and
Kliegl (2003). For \(N=5\):
\[\frac{\sum_{j=1}^{(5-1)/2}x[i+j] -
x[i-j]}{\sum_{j=1}^{(5-1)/2}2j\cdot\Delta t}=\] \[\frac{\sum_{j=1}^{2}x[i+j] -
x[i-j]}{\sum_{j=1}^{2}2j\cdot\Delta t}=\]

\[\frac{x[i+1] - x[i-1]+x[i+2] - x[i-2]}{2\Delta t + 4\Delta t}=\] \[\frac{x[i+2] + x[i+1] - x[i-1] - x[i-2]}{6\Delta t}\]

Horizontal and vertical components of velocity are computed from
adjacent samples, the amplitude is computed from the components and then
is filtered using the Savitzky-Golay
filter. Filter order is controlled via
`nh_sg_filter_order`

option and defaults to
`2`

.

The two methods produce comparable but different velocity values.
However, this has only a minor effect on saccade detection, see a
vignette on *Using sample votes*.

```
# computing velocity using two methods
<- saccadr::diff_ek(x = single_trial$x,
vel_ek y = single_trial$y,
trial = rep(1, nrow(single_trial)),
sample_rate = 500)
<- saccadr::diff_nh(x = single_trial$x,
vel_nh y = single_trial$y,
trial = rep(1, nrow(single_trial)),
sample_rate = 500)
<- na.omit(data.frame(EK = vel_ek[['amp']], NH = vel_nh[['amp']])) %>%
vel_comparison ::filter(NH > 0) %>% # values can be negative due to filtering
dplyr::mutate(logEK = log(EK), logNH = log(NH))
dplyr
# plotting amplitude comparison
<- cor(vel_comparison[['EK']], vel_comparison[['NH']])
pearson_rho
ggplot(data = vel_comparison, aes(x = EK, y = NH)) +
geom_abline() +
geom_point() +
scale_x_log10() +
scale_y_log10() +
labs(subtitle = sprintf("Pearson's ρ = %.2f\nPearson's ρ for log transformed velocities = %.2f", pearson_rho, cor(vel_comparison[['logEK']], vel_comparison[['logNH']])),
x = "Vek",
y = "Vnh")
```

```
# plotting
ggplot(data = NULL, aes(x = vel_comparison[['EK']] - vel_comparison[['NH']])) +
geom_histogram(bins = 50) +
xlab("Vek - Vnh") +
labs(title = "Difference in velocities computed via two methods")
```

A custom function should be as a follows

```
<- function(x, y, trial, sample_rate, options=NULL) {
diff_custom # compute vertical and horizontal components
<- data.frame(x = ..., # differentiate values of x
df y = ...) # differentiate values of y
# compute amplitude
'amp']] = sqrt(df[['x']]^2 + df[['y']]^2))
df[[
# return a data.frame with columns x, y, and amp
return(df)
}
```

Arguments that are passed to the function (all vectors have the same length):

`x`

,`y`

vectors with values, either x and y gaze coordinates or velocity components (for acceleration).`sample_rate`

scalar value in Hz.`trial`

vector with trial index per sample.`options`

a named list with method-specific options. See the example below for how to use them. See also`option_or_default()`

function.

The function must return a data.frame with three columns:
`x`

(horizontal component), `y`

(vertical
component), and `amp`

(amplitude).

Here is an example of a simple velocity function that respects trial borders

```
library(dplyr)
<- function(x, y, trial, sample_rate, options=NULL) {
diff_custom # computing frame time step
<- 1 / sample_rate
delta_t_s
# --- differentiate (compute velocity or acceleration) and filter
data.frame(trial = trial,
x = x,
y = y) %>%
# compute velocity and acceleration for each trial
group_by(trial) %>%
mutate(x = (x - lag(x)) / delta_t_s,
y = (y - lag(y)) / delta_t_s,
amp = sqrt(x^2 + y^2))
}
```

And here is a similar function but that uses an optional
`custom_velocity_lag`

parameter so that it can be different
from 1

```
<- function(x, y, trial, sample_rate, options=NULL) {
diff_lag_n # obtain lag or use default
<- saccadr::option_or_default(options, "custom_velocity_lag", 1)
lag_n
# computing frame time step
<- 1 / sample_rate
delta_t_s
# --- differentiate (compute velocity or acceleration) and filter
data.frame(trial = trial,
x = x,
y = y) %>%
# compute velocity and acceleration for each trial
group_by(trial) %>%
mutate(x = (x - lag(x, lag_n)) / delta_t_s,
y = (y - lag(y, lag_n)) / delta_t_s,
amp = sqrt(x^2 + y^2))
}
```

Engbert, Ralf, and Reinhold Kliegl. 2003. “Microsaccades uncover the orientation of covert
attention.” *Vision Research* 43 (9): 1035–45. https://doi.org/10.1016/S0042-6989(03)00084-1.

Nyström, Marcus, and Kenneth Holmqvist. 2010. “An adaptive algorithm for fixation, saccade, and glissade
detection in eyetracking data.” *Behavior Research
Methods* 42 (1): 188–204. https://doi.org/10.3758/BRM.42.1.188.

Otero-Millan, J., Jose L Alba Castro, Stephen L Macknik, and S.
Martinez-Conde. 2014. “Unsupervised
clustering method to detect microsaccades.” *Journal of
Vision* 14 (2): 18–18. https://doi.org/10.1167/14.2.18.