02-plot-data.Rmd

---
title: "Plot the eyetracking data"
author: "Tristan Mahr"
date: "`r Sys.Date()`"
output: 
  github_document:
    toc: true
    toc_depth: 4
---

```{r setup, include = FALSE}
library("knitr")
opts_chunk$set(
  cache.path = "assets/cache/02-",
  fig.path = "assets/figure/02-",
  warning = FALSE,
  collapse = TRUE, 
  comment = "#>", 
  message = FALSE,
  fig.width = 8,
  fig.asp = 0.618,
  dpi = 300,
  out.width = "80%")
options(width = 100)
```

This script prepares all the various plots we considered using for our
presentation. It also saves the final datasets used for the growth 
curve models.

## Set up

```{r, results = 'hide'}
library(dplyr)
library(littlelisteners)
library(ggplot2)
source("./plotting-helpers.R", encoding = "UTF8")
looks <- readr::read_csv("./data/screened.csv.gz") %>% 
  mutate(
    Cond_Lab = Condition %>% 
      factor(c("real", "MP", "nonsense"),
             c("Real word", "Mispronunciation", "Nonword")))

scores <- readr::read_csv("./data/scores.csv")
```

## Visualize differences in standardized test results

First, confirm that the groups are age-matched. Each step on the _x_-axis
is a pair of children matched by age, sex, and maternal education.

```{r age-match, fig.width = 5, out.width = "60%"}
ggplot(scores) +
  aes(x = resequence(Matching_PairNumber), y = EVT_Age, color = Group_Lab) +
  geom_smooth(method = "lm", formula = y ~ 1, fullrange = TRUE, na.rm = TRUE) +
  # Dodge a bit to avoid drawing points top of each other
  geom_point(na.rm = TRUE, position = position_dodge(.5)) + 
  scale_x_continuous_index() +
  legend_top() + 
  xlab("Pairs matched on age, sex and maternal ed.") +
  ylab("Age (months)") + 
  labs(caption = plot_text$group_intercept_smooth,
       color = "Group")
```

Now, do the same for the standardized test scores. These plots show how the
children systematically differ in vocabulary size and articulation.

```{r scores-by-group, fig.width = 5, out.width = "60%"}
ggplot(scores) +
  aes(x = resequence(Matching_PairNumber), y = EVT_Standard, 
      color = Group_Lab) +
  geom_smooth(method = "lm", formula = y ~ 1, 
              fullrange = TRUE, na.rm = TRUE) +
  geom_point(na.rm = TRUE) + 
  scale_x_continuous_index() +
  legend_top() + 
  align_axis_right() +
  labs(
    x = plot_text$x_match_pairs,
    y = "EVT-2 standard score",
    caption = plot_text$group_intercept_smooth,
    color = "Group")

last_plot() + 
  aes(y = EVT_GSV) + 
  ylab("EVT-2 growth scale value")

last_plot() + 
  aes(y = PPVT_Standard) + 
  ylab("PPVT-4 standard score")

last_plot() + 
  aes(y = GFTA_Standard) + 
  ylab("GFTA-2 standard score")

last_plot() + 
  aes(y = EVT_Age) + 
  ylab("Age (months)")
```

## Visualize the eyetracking data 

### Downsample into 50ms bins

```{r}
times <- looks %>% 
  distinct(Time) %>% 
  arrange(Time) %>% 
  trim_to_bin_width(bin_width = 3, key_time = 0, key_position = 2, Time,
                    min_time = -900) %>% 
  assign_bins(bin_width = 3, Time, bin_col = "Bin") %>% 
  group_by(Bin) %>% 
  mutate(BinTime = median(Time) %>% round(-1))

looks <- looks %>% 
  inner_join(times)
```

### Look at overall looks to familiar image

In order to aggregate looks with my littlelisteners package, we need to define 
a response coding definition so it knows what's a target, what's offscreen, etc.

```{r}
def <- create_response_def(
  label = "LWL scheme",
  primary = "Target",
  others = "Distractor",
  elsewhere = "tracked", 
  missing = NA)
```

We observe that the groups differ on average in the real word and nonword
conditions.

```{r all-conditions-stacked, fig.asp = 1 / 0.618, fig.width = 6}
look_labs <- looks %>% 
  select(Group, Condition, ends_with("_Lab")) %>% 
  distinct()

binned <- looks %>% 
  aggregate_looks(def, Group + Study + ResearchID + 
                    Condition + BinTime ~ GazeByImageAOI) %>% 
  rename(Time = BinTime) %>% 
  left_join(look_labs)

ggplot(binned) + 
  aes(x = Time, y = Prop, color = Group_Lab) + 
  hline_chance() +
  vline_onset() +
  stat_mean_se() + 
  facet_wrap("Cond_Lab", ncol = 1, scales = "free_x") + 
  labs(
    x = plot_text$x_time, 
    y = plot_text$y_target, 
    caption = plot_text$caption_mean_se,
    color = plot_text$group) + 
  legend_top()
```

```{r by-condition}
after_onset <- binned %>% filter(0 <= Time) 

after_onset %>% 
  filter(Group == "NormalHearing") %>% 
  ggplot() + 
    theme_grey(base_size = 12) + 
    aes(x = Time, y = Prop, color = Cond_Lab, shape = Group_Lab) + 
    hline_chance() +
    vline_onset() +
    stat_mean_se(data = after_onset, alpha = 0) + 
    stat_mean_se(data = . %>% filter(Condition == "real")) + 
    viridis::scale_color_viridis(end = .8, discrete = TRUE) + 
    labs(
      x = plot_text$x_time, 
      y = plot_text$y_target, 
      caption = plot_text$caption_mean_se,
      color = "Child hears",
      shape = "Group") + 
    legend_top() + 
    align_axis_right() + 
    expand_limits(x = c(.1, .9))

after_onset %>% 
  ggplot() + 
    theme_grey(base_size = 12) + 
    aes(x = Time, y = Prop, color = Cond_Lab, shape = Group_Lab) + 
    hline_chance() +
    vline_onset() +
    stat_mean_se(data = after_onset, alpha = 0) + 
    stat_mean_se(data = . %>% filter(Condition == "real")) + 
    viridis::scale_color_viridis(end = .8, discrete = TRUE) + 
    labs(
      x = plot_text$x_time, 
      y = plot_text$y_target, 
      caption = plot_text$caption_mean_se,
      color = "Child hears",
      shape = "Group") + 
    legend_top() + 
    align_axis_right() + 
    expand_limits(x = c(.1, .9))

after_onset %>% 
  ggplot() + 
    theme_grey(base_size = 12) + 
    aes(x = Time, y = Prop, color = Cond_Lab, shape = Group_Lab) + 
    hline_chance() +
    vline_onset() +
    stat_mean_se(data = after_onset, alpha = 0) + 
    stat_mean_se(data = . %>% 
                   filter(Condition %in% c("real", "nonsense")) %>% 
                   filter(Group == "NormalHearing")) + 
    stat_mean_se(data = . %>% filter(Condition %in% c("real"))) + 
    viridis::scale_color_viridis(end = .8, discrete = TRUE) + 
    labs(
      x = plot_text$x_time, 
      y = plot_text$y_target, 
      caption = plot_text$caption_mean_se,
      color = "Child hears",
      shape = "Group") + 
    legend_top() + 
    align_axis_right() + 
    expand_limits(x = c(.1, .9))

after_onset %>% 
  ggplot() + 
    theme_grey(base_size = 12) + 
    aes(x = Time, y = Prop, color = Cond_Lab, shape = Group_Lab) + 
    hline_chance() +
    vline_onset() +
    stat_mean_se(data = after_onset, alpha = 0) + 
    stat_mean_se(data = . %>% 
                   filter(Condition %in% c("real", "nonsense")) %>% 
                   filter(Group == "NormalHearing")) + 
    stat_mean_se(data = . %>% filter(Condition %in% c("real", "nonsense"))) + 
    viridis::scale_color_viridis(end = .8, discrete = TRUE) + 
    labs(
      x = plot_text$x_time, 
      y = plot_text$y_target, 
      caption = plot_text$caption_mean_se,
      color = "Child hears",
      shape = "Group") + 
    legend_top() + 
    align_axis_right() + 
    expand_limits(x = c(.1, .9))

after_onset %>% 
  ggplot() + 
    theme_grey(base_size = 12) + 
    aes(x = Time, y = Prop, color = Cond_Lab, shape = Group_Lab) + 
    hline_chance() +
    vline_onset() +
    stat_mean_se(data = after_onset, alpha = 0) + 
    stat_mean_se(data = . %>% 
                   filter(Condition %in% c("real", "nonsense", "MP")) %>% 
                   filter(Group == "NormalHearing")) + 
    stat_mean_se(data = . %>% filter(Condition %in% c("real", "nonsense"))) + 
    viridis::scale_color_viridis(end = .8, discrete = TRUE) + 
    labs(
      x = plot_text$x_time, 
      y = plot_text$y_target, 
      caption = plot_text$caption_mean_se,
      color = "Child hears",
      shape = "Group") + 
    legend_top() + 
    align_axis_right() + 
    expand_limits(x = c(.1, .9))

after_onset %>% 
  ggplot() + 
    theme_grey(base_size = 12) + 
    aes(x = Time, y = Prop, color = Cond_Lab, shape = Group_Lab) + 
    hline_chance() +
    vline_onset() +
    stat_mean_se() + 
    viridis::scale_color_viridis(end = .8, discrete = TRUE) + 
    labs(
      x = plot_text$x_time, 
      y = plot_text$y_target, 
      caption = plot_text$caption_mean_se,
      color = "Child hears",
      shape = "Group") + 
    legend_top() + 
    align_axis_right() + 
    expand_limits(x = c(.1, .9))





after_onset %>% 
  ggplot() + 
    theme_grey(base_size = 12) + 
    aes(x = Time, y = Prop, color = Cond_Lab, shape = Group_Lab) + 
    hline_chance() +
    vline_onset() +
    stat_mean_se(data = after_onset, alpha = 0) + 
    stat_mean_se(data = . %>% 
                   filter(Condition %in% c("nonsense")) %>% 
                   filter(Group == "NormalHearing")) + 
    viridis::scale_color_viridis(end = .8, discrete = TRUE) + 
    labs(
      x = plot_text$x_time, 
      y = plot_text$y_target, 
      caption = plot_text$caption_mean_se,
      color = "Child hears",
      shape = "Group") + 
    legend_top() + 
    align_axis_right() + 
    expand_limits(x = c(.1, .9))

after_onset %>% 
  ggplot() + 
    theme_grey(base_size = 12) + 
    aes(x = Time, y = Prop, color = Cond_Lab, shape = Group_Lab) + 
    hline_chance() +
    vline_onset() +
    stat_mean_se(data = after_onset, alpha = 0) + 
    stat_mean_se(data = . %>% filter(Condition %in% c("nonsense"))) +
    viridis::scale_color_viridis(end = .8, discrete = TRUE) + 
    labs(
      x = plot_text$x_time, 
      y = plot_text$y_target, 
      caption = plot_text$caption_mean_se,
      color = "Child hears",
      shape = "Group") + 
    legend_top() + 
    align_axis_right() + 
    expand_limits(x = c(.1, .9))

after_onset %>% 
  ggplot() + 
    theme_grey(base_size = 12) + 
    aes(x = Time, y = Prop, color = Cond_Lab, shape = Group_Lab) + 
    hline_chance() +
    vline_onset() +
    stat_mean_se(data = after_onset, alpha = 0) + 
    stat_mean_se(data = . %>% 
                   filter(Condition %in% c("MP")) %>% 
                   filter(Group == "NormalHearing")) + 
    viridis::scale_color_viridis(end = .8, discrete = TRUE) + 
    labs(
      x = plot_text$x_time, 
      y = plot_text$y_target, 
      caption = plot_text$caption_mean_se,
      color = "Child hears",
      shape = "Group") + 
    legend_top() + 
    align_axis_right() + 
    expand_limits(x = c(.1, .9))

after_onset %>% 
  ggplot() + 
    theme_grey(base_size = 12) + 
    aes(x = Time, y = Prop, color = Cond_Lab, shape = Group_Lab) + 
    hline_chance() +
    vline_onset() +
    stat_mean_se(data = after_onset, alpha = 0) + 
    stat_mean_se(data = . %>% filter(Condition %in% c("MP"))) +
    viridis::scale_color_viridis(end = .8, discrete = TRUE) + 
    labs(
      x = plot_text$x_time, 
      y = plot_text$y_target, 
      caption = plot_text$caption_mean_se,
      color = "Child hears",
      shape = "Group") + 
    legend_top() + 
    align_axis_right() + 
    expand_limits(x = c(.1, .9))

after_onset %>% 
  ggplot() + 
    theme_grey(base_size = 12) + 
    aes(x = Time, y = Prop, color = Cond_Lab, shape = Group_Lab) + 
    hline_chance() +
    vline_onset() +
    stat_mean_se(data = after_onset, alpha = 0) + 
    stat_mean_se(data = . %>% filter(Condition %in% c("real", "MP"))) +
    viridis::scale_color_viridis(end = .8, discrete = TRUE) + 
    labs(
      x = plot_text$x_time, 
      y = plot_text$y_target, 
      caption = plot_text$caption_mean_se,
      color = "Child hears",
      shape = "Group") + 
    legend_top() + 
    align_axis_right() + 
    expand_limits(x = c(.1, .9))

after_onset %>% 
  ggplot() + 
    theme_grey(base_size = 12) + 
    aes(x = Time, y = Prop, color = Cond_Lab, shape = Group_Lab) + 
    hline_chance() +
    vline_onset() +
    stat_mean_se(data = after_onset, alpha = 0) + 
    stat_mean_se(data = . %>% filter(Condition %in% c("real", "MP"))) +
    viridis::scale_color_viridis(end = .8, discrete = TRUE) + 
    labs(
      x = plot_text$x_time, 
      y = plot_text$y_target, 
      caption = plot_text$caption_mean_se,
      color = "Child hears",
      shape = "Group") + 
    facet_wrap("Group") +
    legend_top() + 
    align_axis_right() + 
    expand_limits(x = c(.1, .9))
```


But given that the groups differ in vocabulary, we can see the same pattern of
results, kind of, by binning children based on EVT scores.

```{r evt-tertile, fig.asp = 1 / 0.618, fig.width = 6}
binned %>% 
  inner_join(scores) %>% 
  filter(!is.na(EVT_Standard)) %>% 
  ggplot() + 
    aes(x = Time, y = Prop, color = factor(ntile(EVT_GSV, 3))) + 
    hline_chance() +
    vline_onset() +
    stat_mean_se() + 
    viridis::scale_color_viridis(discrete = TRUE, end = .7, option = "B") +
    facet_wrap("Cond_Lab", ncol = 1, scales = "free_x") + 
    labs(
      x = plot_text$x_time, 
      y = plot_text$y_target, 
      caption = plot_text$caption_mean_se,
      color = "EVT-2 GSV score tertile") + 
    legend_top() + 
    align_axis_right()
```

The spaghetti plot shows that there are indeed more low-EVT children in the
cochlear implant group.

```{r evt-tertile-spaghetti, fig.asp = 1, fig.width = 6, out.width = "100%"}
binned %>% 
  inner_join(scores) %>% 
  filter(!is.na(EVT_Standard)) %>% 
  filter(0 <= Time) %>% 
  ggplot() + 
    aes(x = Time, y = Prop, color = factor(ntile(EVT_GSV, 3))) + 
    hline_chance() +
    geom_line(aes(group = interaction(Condition, Study, ResearchID))) + 
    viridis::scale_color_viridis(discrete = TRUE, end = .8, option = "B") +
    facet_grid(Cond_Lab ~ Group_Lab, scales = "free_x") + 
    labs(
      x = plot_text$x_time, 
      y = plot_text$y_target, 
      caption = plot_text$caption_mean_se,
      color = "EVT-2 GSV tertile") + 
    legend_top() + 
    align_axis_right()
```

### Get descriptives for slide

```{r}
scores %>% 
  distinct(Group, ChildID, Female) %>% 
  group_by(Group) %>% 
  count(Female)

scores %>% 
  distinct(Group, ChildStudyID) %>% 
  count(Group)

scores %>% 
  group_by(Group) %>% 
  summarise(
    Age = mean(EVT_Age) %>% round(),
    Age_SD = sd(EVT_Age) %>% round(),
    Min_Age = min(EVT_Age),
    Max_Age = max(EVT_Age))

scores %>% 
  group_by(Group) %>% 
  summarise(
    EVT = mean(EVT_Standard, na.rm = TRUE) %>% round(),
    EVT_SD = sd(EVT_Standard, na.rm = TRUE) %>% round(),
    Min_EVT = min(EVT_Standard, na.rm = TRUE),
    Max_EVT = max(EVT_Standard, na.rm = TRUE))

scores %>% 
  distinct(Group, ChildID, Maternal_Education) %>% 
  count(Group, Maternal_Education)
```



### Save the model-ready data

```{r}
zscore <- function(...) as.vector(scale(...))
scores <- scores %>%
  mutate_at(vars(EVT_GSV, EVT_Standard, PPVT_GSV, PPVT_Standard,
         GFTA_Standard), funs(z = zscore))

binned %>% 
  inner_join(scores) %>% 
  select(ChildStudyID, Matching_PairNumber,
         Group:Time, Cond_Lab, Group_Lab, Target, Distractor, Prop,
         EVT_Age, EVT_GSV, EVT_Standard, PPVT_GSV, PPVT_Standard,
         GFTA_Standard, ends_with("_z")) %>% 
  filter(300 <= Time, Time <= 1800) %>% 
  readr::write_csv("./data/model.csv.gz")
```

```{r, include = FALSE}
looks %>% 
  aggregate_looks(def, Bias_ImageAOI + Group + Study + ResearchID + 
                    Condition + BinTime ~ GazeByImageAOI) %>% 
  rename(Time = BinTime) %>% 
  left_join(look_labs) %>% 
  inner_join(scores) %>% 
  select(ChildStudyID, Matching_PairNumber, Bias_ImageAOI,
         Group:Time, Cond_Lab, Group_Lab, Target, Distractor, Prop,
         EVT_Age, EVT_GSV, EVT_Standard, PPVT_GSV, PPVT_Standard,
         GFTA_Standard, ends_with("_z")) %>% 
  filter(300 <= Time, Time <= 1800) %>% 
  readr::write_csv("./data/bias.csv.gz")
```

***

```{r}
sessioninfo::session_info()
```