Using GazeR to Analyze Pupillometry Data

class: center, middle, inverse, title-slide

# Using GazeR to Analyze Pupillometry Data
### Jason Geller
### Rutgers University Center for Cognitive Science

---

# Cognitive Pupillometry

<div class="figure" style="text-align: center">
<img src="amp2020_files/figure-html/unnamed-chunk-2-1.png" alt="PubMed search for the keyword [pupillometry] from 1940-2019"  />
<p class="caption">PubMed search for the keyword [pupillometry] from 1940-2019</p>
</div>
---

# Cognitive Pupillometry

.center[![](assets/img/old.JPG)]

---
# Why Gazer
--

+ GazeR is a package to pre-process data from gaze (e.g., VWP) and pupillometry experiments

+ Not many pipelines available (this is slowly changing)

- PupillometryR (Sam Forbes)
 - Pupillometry (Jason Tsukahara)
 
--

+ The goal was to create a start to finish pipeline for pupillometry

- Work in progress, but create a package that can be used by many different eye tracker brands

+ Amied to simplify some of the pre-processing decisons

+ Much of the psychophysics has been Matlab based, leaving R users out and those who cannot afford Matlab 
---
---
# What does GazeR add?
--

+ Standardized and reproducible pipeline for preprocessing

+ Flexibility in choosing which preprocessing methods and parameters are used

+ Easy to understand pipeline

+ Package can be run in a container through Code Ocean!

- No more worrying about OS or package incompatiablies

---

# Pupil Pre-processing

+ Read in data

+ De-blinking

- Extending blinks
  
+ Interpolation and Smoothing

+ Baseline correction

+ Artifact rejection

- Missing data
  - Unlikely pupil values
  - Median absolute deviation
  
+ Event time alignment

+ Samples to bins (downsampling)

---

# Get Started
--

- Install the **gazeR** package from [Github](https://github.com/dmirman/gazeR):

```r
remotes::install_github("dmirman/gazer")
```

```r
library(gazer)
```

You are recommended to use the [RStudio IDE](https://www.rstudio.com/products/rstudio/), but you do not have to.

---
# Data
+ 4 participants

+ Lexical decison task with easy (type-print) and hard (cursive) trials

+ Data collected with SR EyeLink 1000 plus
---
# Reading in data

+ SR Data Viewer
+ EDF 
+ Tobii (coming soon)

- SR Data Viewer

```r
# set working directory 
setwd("/Users/gellr/Desktop/r")
#folder path 
file_list <- list.files (path="/Users/gellr/Desktop/r", pattern=".xls")
#specify blink, pupil, and whether it is from sr or raw edf files
merge_files=merge_gazer_files(file_list, blink_colname="AVERAGE_IN_BLINK", pupil_colname="AVERAGE_PUPIL_SIZE", filetype="sr")
```
---

# Reading in data

- EDF

```r
directory_edf ="" 
# path to edfs
directory_csv_from_edf_conversion ="" 
# path where csv files should be stored

file_list_edf <- list.files(path=directory_edf, pattern=".edf") 
# get list of edf files
parse_edf(file_list=file_list_edf, 
output_dir = directory_csv_from_edf_conversion, 
type="pupil")) 
# parses edfs and saves them into directory
# folder path to csv folders from parse_edf

file_list_pupil_samp <- list.files(path=directory_csv_from_edf_conversion, 
pattern=".csv") 
# extract the processed csv files from directory_csv_from_edf_conversion

pd <- merge_gazer_files(file_list_pupil_samp, type = "edf") 
# merges all the files from file_list_pupil_samp object
```
---
# Reading in Data Example

```r
pupil_path <- system.file("extdata", "pupil_sample_files_edf.xls", package = "gazer")
pupil_raw<-fread(pupil_path)
pupil_raw<-as_tibble(pupil_raw)

knitr::kable(head(pupil_raw))
```

|subject | trial| time| pupil|     x|     y| blink|message       | acc| block|   rt|item      |script  |alt  |
|:-------|-----:|----:|-----:|-----:|-----:|-----:|:-------------|---:|-----:|----:|:---------|:-------|:----|
|11c.edf |     1|    0|  3635| 971.6| 563.2|     0|MODERECORDCRL |   1|     0| 3833|mourn.png |cursive |word |
|11c.edf |     1|    4|  3634| 973.4| 564.7|     0|NA            |   1|     0| 3833|mourn.png |cursive |word |
|11c.edf |     1|    8|  3626| 971.8| 566.7|     0|NA            |   1|     0| 3833|mourn.png |cursive |word |
|11c.edf |     1|   12|  3622| 971.0| 566.2|     0|NA            |   1|     0| 3833|mourn.png |cursive |word |
|11c.edf |     1|   16|  3621| 972.1| 566.2|     0|NA            |   1|     0| 3833|mourn.png |cursive |word |
|11c.edf |     1|   20|  3621| 972.1| 565.8|     0|NA            |   1|     0| 3833|mourn.png |cursive |word |

---
# Data Structure

+ Can work with with data from any tracker as long as you have these columns:

```r
knitr::kable(names_needed, col.names = "Names")
```

|Names      |
|:----------|
|Subject    |
|Trial      |
|Blink      |
|Pupil Size |
|Time       |
|Message    |
---
# De-Blinking

.top[![](assets/img/image5.png)]

---

# De-Blinking

+ Blinks are a major artifact in pupillometry research

- Eye-trackers treat blinks as missing data
 
+ Values before and after blinks are spurious

+ We need to extend them
--

```r
pup_extend<- pupil_raw %>% 
  group_by(subject, trial) %>% 
  mutate(extendpupil=extend_blinks(pupil, fillback=100, 
                                   fillforward=100, hz=250))
  
pup_extend1 <- pup_extend %>%
  select(subject, trial, time, pupil, extendpupil)
```

---

# Smoothing/Filtering and Interpolation

+ Do not throw out data!

+ Let's interpolate and smooth the data instead

- GazeR can do linear or cubic interpolation 
--

+ Smoothing with n moving average or Hanning filter

+ Choice seems to matter
   
--

```r
smooth_interp <- smooth_interpolate_pupil(pup_extend, pupil="pupil", 
                                          extendpupil="extendpupil", 
                                          extendblinks=TRUE, step.first="interp",                                                                      filter="moving", 
                                          maxgap=Inf, type="linear", hz=250, n=5)
```

```
## Performing linear interpolation
```

```
## Smoothing the pupil trace with moving average
```

---

<img src="assets/img/lin.png" width="65%" style="display: block; margin: auto;" /><img src="assets/img/cubic.png" width="65%" style="display: block; margin: auto;" />
---
# Baseline Correction

+ Subtractive or divisive (grand average baseline)

```r
baseline_pupil <- baseline_correction_pupil(smooth_interp, 
                                            pupil_colname='pup_interp', 
                                            baseline_window=c(500,1000))

baseline_pupil<-baseline_correction_pupil_msg(smooth_interp, pupil_colname='pup_interp', 
                                              baseline_dur=100, event="target", 
                                              baseline_method = "sub")

knitr::kable(head(baseline_pupil))
```

|subject | trial| baseline| time| pupil|     x|     y| blink|message       | acc| block|   rt|item      |script  |alt  | extendpupil| interp| pup_interp| baselinecorrectedp|
|:-------|-----:|--------:|----:|-----:|-----:|-----:|-----:|:-------------|---:|-----:|----:|:---------|:-------|:----|-----------:|------:|----------:|------------------:|
|11c.edf |     1|   3584.3|    0|  3635| 971.6| 563.2|     0|MODERECORDCRL |   1|     0| 3833|mourn.png |cursive |word |        3635|   3635|   3631.667|           47.36667|
|11c.edf |     1|   3584.3|    4|  3634| 973.4| 564.7|     0|NA            |   1|     0| 3833|mourn.png |cursive |word |        3634|   3634|   3629.250|           44.95000|
|11c.edf |     1|   3584.3|    8|  3626| 971.8| 566.7|     0|NA            |   1|     0| 3833|mourn.png |cursive |word |        3626|   3626|   3627.600|           43.30000|
|11c.edf |     1|   3584.3|   12|  3622| 971.0| 566.2|     0|NA            |   1|     0| 3833|mourn.png |cursive |word |        3622|   3622|   3624.800|           40.50000|
|11c.edf |     1|   3584.3|   16|  3621| 972.1| 566.2|     0|NA            |   1|     0| 3833|mourn.png |cursive |word |        3621|   3621|   3621.600|           37.30000|
|11c.edf |     1|   3584.3|   20|  3621| 972.1| 565.8|     0|NA            |   1|     0| 3833|mourn.png |cursive |word |        3621|   3621|   3620.000|           35.70000|
---

# Artifact Rejection: Missingness

+ Should remove trials and subjects with more than 20-30% of data

```r
pup_missing <- count_missing_pupil(baseline_pupil, pupil="pupil", 
                                   missingthresh = .2)
```

```
## % trials excluded:0.079
```

```
## Participants taken out:character(0)
```

---
# Artifact Rejection: Spurious Values

+ Avoid using heuristics
  - +/- SD
  
+ Visualize your data

+ Small baseline values can exaggerate differences (percent change)

```r
puphist <- ggplot(baseline_pupil, aes(x = pup_interp)) + 
  geom_histogram(aes(y = ..count..), colour = "green", binwidth = 0.5)  + 
  geom_vline(xintercept = 2500, linetype="dotted") +
  geom_vline(xintercept = 4000, linetype="dotted") + 
  xlab("Pupil Size") + 
  ylab("Count") + 
  theme_bw()
```
---
# Histogram

```r
puphist
```

```r
pup_outliers <- pup_missing %>% 
  # based on visual inspection
  dplyr::filter(pup_interp  >= 2500, pup_interp <= 4000) 
```

---

# Artifact rejection: MAD

$$
{d}^{\prime \left[i\right]}=\max \left(\left|\frac{d\left[i\right]-d\left[i-1\right]}{t\left[i\right]-t\left[i-1\right]}\right|,\left|\frac{d\left[i+1\right]-d\left[i\right]}{t\left[i+1\right]-t\left[i\right]}\right|\right)
$$

$$
\mathrm{MAD}=\mathrm{median}\left(\left|{d}^{\prime}\left[i\right]-\mathrm{median}\left({d}^{\prime}\right)\right|\right)
$$

$$
\mathrm{Treshold}=\mathrm{median}\left({d}^{\prime}\right)+n\bullet \mathrm{MAD}
$$

```r
mad_removal <- pup_outliers  %>% 
  group_by(subject, trial) %>% 
  mutate(speed=speed_pupil(pup_interp,time)) %>% 
  mutate(MAD=calc_mad(speed, n = 16)) %>% 
  filter(speed < MAD)
```

---
# Event time Alignment

```r
baseline_pupil_onset <- mad_removal %>% 
  dplyr::group_by(subject, trial) %>%  
  dplyr::mutate(time_zero=onset_pupil(time, message, event=c("target"))) %>%
  ungroup() %>% 
  dplyr::filter(time_zero >= 0 & time_zero <= 2500) %>%
  select(subject, trial, time_zero, message, baselinecorrectedp, pup_interp, script)

knitr::kable(head(baseline_pupil_onset))
```

|subject | trial| time_zero|message | baselinecorrectedp| pup_interp|script  |
|:-------|-----:|---------:|:-------|------------------:|----------:|:-------|
|11c.edf |     1|         0|target  |               -4.1|     3580.2|cursive |
|11c.edf |     1|         4|NA      |               -4.1|     3580.2|cursive |
|11c.edf |     1|         8|NA      |               -3.9|     3580.4|cursive |
|11c.edf |     1|        12|NA      |               -4.5|     3579.8|cursive |
|11c.edf |     1|        16|NA      |               -4.7|     3579.6|cursive |
|11c.edf |     1|        20|NA      |               -5.3|     3579.0|cursive |

---
# Samples to bins

```r
timebins1<- gazer::downsample_gaze(baseline_pupil_onset, bin.length=100, timevar = "time_zero", aggvars = c("subject", "script", "timebins"), type="pupil")
```

```
## Warning: The `.dots` argument of `group_by()` is deprecated as of dplyr 1.0.0.
```

```r
knitr::kable(head(timebins1))
```

|subject |script  | timebins| aggbaseline|
|:-------|:-------|--------:|-----------:|
|11c.edf |cursive |        0|   -2.869280|
|11c.edf |cursive |      100|   -3.555438|
|11c.edf |cursive |      200|   -4.038295|
|11c.edf |cursive |      300|   -1.726396|
|11c.edf |cursive |      400|    2.292393|
|11c.edf |cursive |      500|    4.686093|

---
# Visualize

---

# Analysis: GCA

+ gazer includes a helper fucntion to perfrom GCA

```r
code_poly(df = timebins1, predictor ="timebins", poly.order =3, orthogonal =TRUE, draw.poly =TRUE)
```

```
## # A tibble: 156 x 8
##    subject script  timebins aggbaseline timebins.Index   poly1     poly2   poly3
##    <chr>   <chr>      <dbl>       <dbl>          <dbl>   <dbl>     <dbl>   <dbl>
##  1 11c.edf cursive        0       -2.87              1 -0.327   3.91e- 1 -0.412 
##  2 11c.edf cursive      100       -3.56              2 -0.301   2.97e- 1 -0.214 
##  3 11c.edf cursive      200       -4.04              3 -0.275   2.11e- 1 -0.0576
##  4 11c.edf cursive      300       -1.73              4 -0.248   1.33e- 1  0.0612
##  5 11c.edf cursive      400        2.29              5 -0.222   6.25e- 2  0.146 
##  6 11c.edf cursive      500        4.69              6 -0.196  -6.82e-17  0.200 
##  7 11c.edf cursive      600       -2.86              7 -0.170  -5.47e- 2  0.228 
##  8 11c.edf cursive      700      -10.1               8 -0.144  -1.02e- 1  0.232 
##  9 11c.edf cursive      800       -8.49              9 -0.118  -1.41e- 1  0.217 
## 10 11c.edf cursive      900       -3.05             10 -0.0915 -1.72e- 1  0.185 
## # … with 146 more rows
```
---

# The Future of gazeR

+ Submit to CRAN

+ Support for more ETs

+ GUI/Shiny

+ Better error and warning messages

+ What do folks want to see?

---

# Thanks!

---