03-models.Rmd

---
title: "Growth curve models"
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/03-",
  fig.path = "assets/figure/03-",
  warning = FALSE,
  collapse = TRUE, 
  comment = "#>", 
  message = FALSE,
  fig.width = 8,
  fig.asp = 0.618,
  dpi = 300,
  out.width = "80%")
options(width = 100)
```

This script fits the growth curve models for the analyses.

## Set up

```{r, results = 'hide'}
library(dplyr)
library(littlelisteners)
library(ggplot2)
library(lme4)

source("./plotting-helpers.R", encoding = "UTF8")
looks <- readr::read_csv("./data/model.csv.gz") %>% 
  mutate(
    Cond_Lab = Condition %>% 
      factor(c("real", "MP", "nonsense"),
             c("Real word", "Mispronunciation", "Nonword")))
```

### Add orthogonal polynomials

```{r}
looks <- looks %>% 
  polypoly::poly_add_columns(Time, degree = 3, prefix = "ot") 
```

### Prep the datasets

For the models with vocabulary information, we need to remove any pairs of
children in which one of the children is missing an EVT score. We identify those 
pairs.

```{r}
no_vocab_pairs <- looks %>%
  distinct(ChildStudyID, Matching_PairNumber, EVT_Standard) %>%
  filter(is.na(EVT_Standard)) %>%
  select(Matching_PairNumber) %>%
  print()
```

Prepare a dataset for each condition.

```{r}
d_mp <- looks %>% 
  filter(Condition == "MP")

d_rw <- looks %>% 
  filter(Condition == "real")

d_ns <- looks %>% 
  filter(Condition == "nonsense")

d_mp_evt <- d_mp %>% 
  anti_join(no_vocab_pairs, by = "Matching_PairNumber")

d_rw_evt <- d_rw %>% 
  anti_join(no_vocab_pairs, by = "Matching_PairNumber")

d_ns_evt <- d_ns %>% 
  anti_join(no_vocab_pairs, by = "Matching_PairNumber")
```




## Fit the vocabulary-less models

First, we just look at the group effects on the growth curve shapes.

```{r no-evt-models, cache = TRUE}
glmer_controls <- glmerControl(
  optimizer = "bobyqa",
  optCtrl = list(maxfun = 2e5))

m_mp <- glmer(
  cbind(Target, Distractor) ~ 
    Group * (ot1 + ot2 + ot3) + (ot1 + ot2 + ot3 | ChildStudyID),
  family = binomial,
  control = glmer_controls,
  data = d_mp)
summary(m_mp)

m_rw <- glmer(
  cbind(Target, Distractor) ~ 
    Group * (ot1 + ot2 + ot3) + (ot1 + ot2 + ot3 | ChildStudyID),
  family = binomial,
  control = glmer_controls,
  data = d_rw)
summary(m_rw)

m_ns <- glmer(
  cbind(Target, Distractor) ~ 
    Group * (ot1 + ot2 + ot3) + (ot1 + ot2 + ot3 | ChildStudyID),
  family = binomial,
  control = glmer_controls,
  data = d_ns)
summary(m_ns)
```

We can plot the growth curve fixed effects which describe how the average child
in each group x condition performs.

```{r overall-fits, echo = FALSE}
d_times <- d_ns %>% 
  distinct(Time, ot1, ot2, ot3, Group)

condition_labels <- looks %>% 
  distinct(Condition, Cond_Lab, Group, Group_Lab)

models <- list(
  nonsense = m_ns,
  MP = m_mp,
  real = m_rw)

fitted <- map_predict(d_times, models, predict_fixef_resp) %>% 
  left_join(condition_labels, 
            by = c("model" = "Condition", "Group"))

ggplot(looks) +
  aes(x = Time, color = Cond_Lab, group = interaction(Group_Lab, Cond_Lab)) + 
  hline_chance() +
  stat_mean_se(mapping = aes(y = Prop), alpha = .4) + 
  geom_line(aes(y = fitted, linetype = Group_Lab, 
                group = interaction(Group_Lab, Cond_Lab)), data = fitted, 
            size = 1) +
  viridis::scale_color_viridis(discrete = TRUE, option = "viridis", end = .8) +
  guides(colour = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = plot_text$x_time,
       y = plot_text$y_fits,
       color = "Child hears",
       linetype = "Group") + 
  legend_top() + 
  align_axis_right()

ggplot(fitted) + 
  aes(x = Time, y = fitted, linetype = Group_Lab, color = Cond_Lab) + 
  hline_chance() +
  geom_line(size = 1) +
  labs(x = plot_text$x_time, y = plot_text$y_fits,
       color = "Child hears", linetype = "Group") + 
  legend_top() + 
  align_axis_right() + 
  viridis::scale_color_viridis(discrete = TRUE, option = "viridis", end = .8)

# logit_fits <- map_predict(d_times, models, predict_fixef_link) %>% 
#   left_join(condition_labels, 
#             by = c("model" = "Condition", "Group"))
# 
# ggplot(logit_fits) + 
#   aes(x = Time, y = fitted, linetype = Group_Lab, color = Cond_Lab) + 
#   hline_chance() +
#   geom_line(size = 1) +
#   labs(x = plot_text$x_time, y = plot_text$y_fits,
#        color = "Child hears", linetype = "Group") + 
#   legend_top() + 
#   align_axis_right()
```




## Fit the models that include vocabulary isze


### Mispronunciations

<!-- The two groups significantly differ with respect to their intercept terms, and -->
<!-- not the shape of their growth curves. -->

<!-- Now, we have to refit the model the model to exclude _pairs_ of children where one of the children is missing an EVT2 score vocabulary scores. -->

Allow EVT to interact with time, time^2^ and time^3^ because these data are 
curvier.

```{r mp-evt-models, cache = TRUE}
m_mp_1a <- glmer(
  cbind(Target, Distractor) ~ 
    Group * (1 + ot1 + ot2 + ot3) + (ot1 + ot2 + ot3 | ChildStudyID),
  family = binomial,
  control = glmer_controls,
  data = d_mp_evt)
summary(m_mp_1a)

m_mp_2a <- glmer(
  cbind(Target, Distractor) ~ 
    EVT_GSV_z * (1 + ot1 + ot2 + ot3) + (ot1 + ot2 + ot3 | ChildStudyID),
  family = binomial,
  control = glmer_controls,
  data = d_mp_evt)
summary(m_mp_2a)

m_mp_3a <- glmer(
  cbind(Target, Distractor) ~ 
    Group * (1 + ot1 + ot2 + ot3) + 
    EVT_GSV_z * (1 + ot1 + ot2 + ot3) + 
    (ot1 + ot2 + ot3 | ChildStudyID),
  family = binomial,
  control = glmer_controls,
  data = d_mp_evt)
summary(m_mp_3a)

m_mp_3b <- glmer(
  cbind(Target, Distractor) ~ 
    Group * EVT_GSV_z * (1 + ot1 + ot2 + ot3) + 
    (ot1 + ot2 + ot3 | ChildStudyID),
  family = binomial,
  control = glmer_controls,
  data = d_mp_evt)
summary(m_mp_3b)
```

The last two plots are from the Group + Vocab models.

```{r mp-plots, fig.width = 6, echo = FALSE}
scale_plasma <- viridis::scale_color_viridis(
  discrete = TRUE, 
  option = "plasma", 
  end = .7)

models <- list(
  Group = m_mp_1a, 
  `Vocab.` = m_mp_2a, 
  `Group + Vocab.` = m_mp_3a, 
  `Group x Vocab.` = m_mp_3b)

d_times <- d_ns %>% 
  distinct(Time, ot1, ot2, ot3, Group)

evt_lmh <- d_mp_evt %>% 
  distinct(Group, Group_Lab, EVT_GSV_z) %>% 
  group_by(Group, Group_Lab) %>% 
  tidy_quantile(EVT_GSV_z)

fits <- d_times %>% 
  inner_join(evt_lmh, by = "Group") %>% 
  map_predict(models, predict_fixef_resp) %>% 
  rename(Model = model) %>% 
  mutate(Model = factor(Model, c("Group", "Vocab.", 
                                 "Group + Vocab.", "Group x Vocab.")))

ggplot(fits) +
  aes(x = Time, y = fitted, color = Group_Lab) + 
  hline_chance() +
  geom_line(
    aes(x = Time, y = fitted, group = interaction(Group, quantile))) + 
  scale_plasma +
  facet_wrap("Model", labeller = label_both) + 
  labs(caption = "Lines: Within group EVT quantiles (10%, 30%, 50%, 70%, 90%)",
     x = plot_text$x_time,
     y = plot_text$y_fits,
     color = "Group") + 
  legend_top() + 
  align_axis_right() +
  expand_limits(y = c(.25, .75))
  
ggplot(fits %>% filter(Model == "Group + Vocab.")) +
  aes(x = Time, y = fitted) + 
  hline_chance() +
  geom_line(
    aes(x = Time, y = fitted, group = interaction(Group, quantile)),
    size = 1) + 
  facet_wrap("Group_Lab") + 
  labs(caption = "Lines: Within group EVT quantiles (10%, 30%, 50%, 70%, 90%)",
       x = plot_text$x_time,
       y = plot_text$y_fits) + 
  legend_top() + 
  align_axis_right() + 
  ggtitle("Vocabulary effects in mispronunciation condition") +
  expand_limits(y = c(.25, .75))

ggplot(fits %>% filter(Model == "Group + Vocab.")) +
  aes(x = Time, y = fitted, color = quantile) + 
  hline_chance() +
  geom_line(
    aes(x = Time, y = fitted, group = interaction(Group, quantile)),
    size = 1) + 
  facet_wrap("Group_Lab") + 
  viridis::scale_color_viridis(discrete = TRUE, end = .85, option = "B") +
  labs(x = plot_text$x_time,
       y = plot_text$y_fits,
       color = "Within group EVT quantiles") + 
  legend_bottom(legend.justification = "right") + 
  align_axis_right() + 
  ggtitle("Vocabulary effects in mispronunciation condition") +
  expand_limits(y = c(.25, .75))

d_mp %>% 
  anti_join(no_vocab_pairs) %>% 
  left_join(condition_labels) %>% 
  group_by(Group) %>% 
  mutate(EVT_Bin = ntile(EVT_GSV, 5) %>% factor()) %>% 
  ungroup() %>% 
  ggplot() +
    aes(x = Time, y = Prop, color = EVT_Bin) + 
    hline_chance() +
    stat_mean_se() +
    facet_wrap("Group_Lab") + 
    viridis::scale_color_viridis(discrete = TRUE, end = .85, option = "B") +
    labs(x = plot_text$x_time,
         y = plot_text$y_fits,
         color = "Within group EVT bins") + 
    legend_bottom(legend.justification = "right") + 
    align_axis_right() + 
    ggtitle("Vocabulary effects in mispronunciation condition") +
  expand_limits(y = c(.25, .75))
```

Model comparison here should be taken with a grain of salt because we are adding
predictors in batches of four coefficients (the predictor's main effect and
interactions with time). AIC and BIC penalize fit measures based on number of
parameters, so we are adding big penalities with each additional predictor.
Consider this more of coarse look at how batches of coefficients improve model
fit.

```{r mp-anova}
do.call(anova, unname(models))
```


### Real words

Allow EVT to interact with time^1^, time^2^ and time^3^.

```{r rw-evt-models, cache = TRUE, eval = TRUE}
m_rw_1a <- glmer(
  cbind(Target, Distractor) ~ 
    Group * (1 + ot1 + ot2 + ot3) + (ot1 + ot2 + ot3 | ChildStudyID),
  family = binomial,
  control = glmer_controls,
  data = d_rw %>% anti_join(no_vocab_pairs))
summary(m_rw_1a)

m_rw_2a <- glmer(
  cbind(Target, Distractor) ~ 
   EVT_GSV_z * (1 + ot1 + ot2 + ot3) + ot2 + ot3 + (ot1 + ot2 + ot3 | ChildStudyID),
  family = binomial,
  control = glmer_controls,
  data = d_rw %>% anti_join(no_vocab_pairs))
summary(m_rw_2a)

m_rw_3a <- glmer(
  cbind(Target, Distractor) ~ 
    Group * (1 + ot1 + ot2 + ot3) + 
    EVT_GSV_z * (1 + ot1 + ot2 + ot3) + 
    (ot1 + ot2 + ot3 | ChildStudyID),
  family = binomial,
  control = glmer_controls,
  data = d_rw %>% anti_join(no_vocab_pairs))
summary(m_rw_3a)

m_rw_3b <- glmer(
  cbind(Target, Distractor) ~ 
    Group * EVT_GSV_z * (1 + ot1 + ot2 + ot3) + 
    (ot1 + ot2 + ot3 | ChildStudyID),
  family = binomial,
  control = glmer_controls,
  data = d_rw %>% anti_join(no_vocab_pairs))
summary(m_rw_3b)
```


The last two plots are from the Group + Vocab model.

```{r rw-plots, fig.width = 6, echo = FALSE}
models <- list(
  Group = m_rw_1a, 
  `Vocab.` = m_rw_2a, 
  `Group + Vocab.` = m_rw_3a, 
  `Group x Vocab.` = m_rw_3b)

fits <- d_times %>% 
  inner_join(evt_lmh, by = "Group") %>% 
  map_predict(models, predict_fixef_resp) %>% 
  rename(Model = model) %>% 
  mutate(Model = factor(Model, c("Group", "Vocab.", 
                                 "Group + Vocab.", "Group x Vocab.")))

ggplot(fits) +
  aes(x = Time, y = fitted, color = Group_Lab) + 
  hline_chance() +
  geom_line(
    aes(x = Time, y = fitted, group = interaction(Group, quantile))) + 
  scale_plasma + 
  facet_wrap("Model", labeller = label_both) + 
    labs(caption = "Lines: Within group EVT quantiles (10%, 30%, 50%, 70%, 90%)",
       x = plot_text$x_time,
       y = plot_text$y_fits,
       color = "Group") + 
  legend_top() + 
  align_axis_right() +
  expand_limits(y = c(.25, .75))
  
ggplot(fits %>% filter(Model == "Group + Vocab.")) +
  aes(x = Time, y = fitted) + 
  hline_chance() +
  geom_line(
    aes(x = Time, y = fitted, group = interaction(Group, quantile)),
    size = 1) + 
  facet_wrap("Group_Lab") + 
  labs(caption = "Lines: Within group EVT quantiles (10%, 30%, 50%, 70%, 90%)",
       x = plot_text$x_time,
       y = plot_text$y_fits) + 
  legend_top() + 
  align_axis_right() + 
  ggtitle("Vocabulary effects in correct production condition")

ggplot(fits %>% filter(Model == "Group + Vocab.")) +
  aes(x = Time, y = fitted, color = quantile) + 
  hline_chance() +
  geom_line(
    aes(x = Time, y = fitted, group = interaction(Group, quantile)),
    size = 1) + 
  facet_wrap("Group_Lab") + 
  viridis::scale_color_viridis(discrete = TRUE, end = .85, option = "B") +
  labs(x = plot_text$x_time,
       y = plot_text$y_fits,
       color = "Within group EVT quantiles") + 
  legend_bottom(legend.justification = "right") + 
  align_axis_right() + 
  ggtitle("Vocabulary effects in correct production condition")
```

This final plot is from the Group x Vocab model. I think the Normal Hearing side
is overly sensitive to some outlying behavior in the high EVT normal hearing
group.

```{r rw-plots-2}
ggplot(fits %>% filter(Model == "Group x Vocab.")) +
  aes(x = Time, y = fitted, color = quantile) + 
  hline_chance() +
  geom_line(
    aes(x = Time, y = fitted, group = interaction(Group, quantile)),
    size = 1) + 
  facet_wrap("Group_Lab") + 
  viridis::scale_color_viridis(discrete = TRUE, end = .85, option = "B") +
  labs(x = plot_text$x_time,
       y = plot_text$y_fits,
       color = "Within group EVT quantiles") + 
  legend_bottom(legend.justification = "right") + 
  align_axis_right() + 
  ggtitle("Vocabulary effects in correct production condition")

d_rw %>% 
  anti_join(no_vocab_pairs) %>% 
  left_join(condition_labels) %>% 
  group_by(Group) %>% 
  mutate(EVT_Bin = ntile(EVT_GSV, 5) %>% factor()) %>% 
  ungroup() %>% 
  ggplot() +
    aes(x = Time, y = Prop, color = EVT_Bin) + 
    hline_chance() +
    stat_mean_se() +
    facet_wrap("Group_Lab") + 
    viridis::scale_color_viridis(discrete = TRUE, end = .85, option = "B") +
    labs(x = plot_text$x_time,
         y = plot_text$y_fits,
         color = "Within group EVT bins") + 
    legend_bottom(legend.justification = "right") + 
    align_axis_right() + 
    ggtitle("Vocabulary effects in correct production condition")
```


```{r real-word-anova}
do.call(anova, unname(models))
```


### Nonword models

Allow EVT to interact with time^1^.

```{r ns-evt-models, cache = TRUE}
m_ns_1a <- glmer(
  cbind(Target, Distractor) ~ 
    Group * (1 + ot1) + ot2 + ot3 + (ot1 + ot2 + ot3 | ChildStudyID),
  family = binomial,
  control = glmer_controls,
  data = d_ns %>% anti_join(no_vocab_pairs))
summary(m_ns_1a)

m_ns_2a <- glmer(
  cbind(Target, Distractor) ~ 
    EVT_GSV_z * (1 + ot1) + ot2 + ot3 + (ot1 + ot2 + ot3 | ChildStudyID),
  family = binomial,
  control = glmer_controls,
  data = d_ns %>% anti_join(no_vocab_pairs))
summary(m_ns_2a)

m_ns_3a <- glmer(
  cbind(Target, Distractor) ~ 
    Group * (1 + ot1 + ot2 + ot3) + 
    EVT_GSV_z * (1 + ot1) + 
    (ot1 + ot2 + ot3 | ChildStudyID),
  family = binomial,
  control = glmer_controls,
  data = d_ns %>% anti_join(no_vocab_pairs))
summary(m_ns_3a)

m_ns_3b <- glmer(
  cbind(Target, Distractor) ~ 
    Group * (1 + ot1 + ot2 + ot3) + 
    Group * EVT_GSV_z * (1 + ot1) + 
    (ot1 + ot2 + ot3 | ChildStudyID),
  family = binomial,
  control = glmer_controls,
  data = d_ns %>% anti_join(no_vocab_pairs))
summary(m_ns_3b)
```

The last two plots are from the Group + Vocab model.

```{r ns-plots, fig.width = 6, echo = FALSE}
models <- list(
  Group = m_ns_1a, 
  `Vocab.` = m_ns_2a, 
  `Group + Vocab.` = m_ns_3a, 
  `Group x Vocab.` = m_ns_3b)

fits <- d_times %>% 
  inner_join(evt_lmh, by = "Group") %>% 
  map_predict(models, predict_fixef_resp) %>% 
  rename(Model = model) %>% 
  mutate(Model = factor(Model, c("Group", "Vocab.", 
                                 "Group + Vocab.", "Group x Vocab.")))
ggplot(fits) +
  aes(x = Time, y = fitted, color = Group_Lab) + 
  hline_chance() +
  geom_line(
    aes(x = Time, y = fitted, group = interaction(Group, quantile))) + 
  scale_plasma + 
  facet_wrap("Model", labeller = label_both) + 
    labs(caption = "Lines: Within group EVT quantiles (10%, 30%, 50%, 70%, 90%)",
       x = plot_text$x_time,
       y = plot_text$y_fits,
       color = "Group") + 
  legend_top() + 
  align_axis_right() +
  expand_limits(y = c(.25, .75))
  
ggplot(fits %>% filter(Model == "Group + Vocab.")) +
  aes(x = Time, y = fitted) + 
  hline_chance() +
  geom_line(
    aes(x = Time, y = fitted, group = interaction(Group_Lab, quantile)),
    size = 1) + 
  facet_wrap("Group_Lab") + 
  labs(caption = "Lines: Within group EVT quantiles (10%, 30%, 50%, 70%, 90%)",
       x = plot_text$x_time,
       y = plot_text$y_fits) + 
  legend_top() + 
  align_axis_right() + 
  ggtitle("Vocabulary effects in nonword condition")

ggplot(fits %>% filter(Model == "Group + Vocab.")) +
  aes(x = Time, y = fitted, color = quantile) + 
  hline_chance() +
  geom_line(
    aes(x = Time, y = fitted, group = interaction(Group, quantile)),
    size = 1) + 
  facet_wrap("Group_Lab") + 
  viridis::scale_color_viridis(discrete = TRUE, end = .85, option = "B") +
  labs(x = plot_text$x_time,
       y = plot_text$y_fits,
       color = "Within group EVT quantiles") + 
  legend_bottom(legend.justification = "right") + 
  align_axis_right() + 
  ggtitle("Vocabulary effects in nonword condition")

d_ns %>% 
  anti_join(no_vocab_pairs) %>% 
  left_join(condition_labels) %>% 
  group_by(Group) %>% 
  mutate(EVT_Bin = ntile(EVT_GSV, 5) %>% factor()) %>% 
  ungroup() %>% 
  ggplot() +
    aes(x = Time, y = Prop, color = EVT_Bin) + 
    hline_chance() +
    stat_mean_se() +
    facet_wrap("Group_Lab") + 
    viridis::scale_color_viridis(discrete = TRUE, end = .85, option = "B") +
    labs(x = plot_text$x_time,
         y = plot_text$y_fits,
         color = "Within group EVT bins") + 
    legend_bottom(legend.justification = "right") + 
    align_axis_right() + 
    ggtitle("Vocabulary effects in nonword condition")
```


```{r nonword-anova}
do.call(anova, unname(models))
```


***


```{r vocab-boxplot, fig.width = 5}
d_ns %>% 
  anti_join(no_vocab_pairs) %>% 
  distinct(Group_Lab, ResearchID, ChildStudyID, EVT_GSV) %>% 
  ggplot() + 
    theme(axis.text.x = element_text(size = rel(1))) + 
    aes(x = Group_Lab, y = EVT_GSV) + 
    geom_boxplot(width = .4) +
    geom_point(position = position_jitter(.1), shape = 1) + 
    labs(x = NULL, y = "EVT-2 GSV") + 
    ggtitle("Children with CIs have more variable vocabularies") + 
    align_axis_right()
```


***

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