04-gca-models.Rmd

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

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

wd <- rprojroot::find_rstudio_root_file()
opts_knit$set(root.dir = wd)
```

## Setup

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

# Work relative to RStudio project
wd <- rprojroot::find_rstudio_root_file()
d_m <- readr::read_csv(file.path(wd, "data", "modeling.csv"))
```

```{r, echo = FALSE}
# These are plotting settings and helpers.

stat_mean_se <- function(...) {
  stat_summary(fun.data = mean_se, ...)
}

stat_mean <- function(...) {
  dots <- list(...)
  dots$fun.y <- mean
  if (is.null(dots$geom)) dots$geom <- "line"
  do.call(stat_summary, dots)
}
  
plot_text <- list(
  x_time = "Time (ms) after target noun onset", 
  y_target = "Proportion of looks to named image",
  y_image = "Proportion of looks to image",
  caption_mean_se = "Mean ± SE"
)

legend_bottom <- theme(
  legend.position = "bottom", 
  legend.text = element_text(size = 10), 
  legend.justification = "left")

legend_top <- theme(
  legend.position = "top", 
  legend.text = element_text(size = 10), 
  legend.justification = "left")

colors <- viridis::scale_color_viridis(end = .7, discrete = TRUE)

hline_chance <- geom_hline(yintercept = .25, size = 1.25, color = "#cccccc")
vline_onset <- geom_vline(xintercept = 0, size = 1.25, color = "#cccccc")
```



## Growth curve models prep

Based on the [visual exploratory data analysis](./03-eyetracking-data.md), there
is no average effect of hearing a non-native dialect on familiar word
recognition. We need to confirm this observation through growth curve analysis.


### Set analysis window

First, we set some options for the analysis, namely the start and end time of
the analysis window.

```{r}
opts_gca <- list(
  min_time = 250,
  max_time = 1500
)
```

We model looks from `r opts_gca$min_time` to `r opts_gca$max_time` ms.

```{r windows, echo = FALSE}
# Split data into three dataframes, based on whether looks are before, 
# after, or inside analysis window
data_windows <- d_m %>% 
  mutate(
    window = case_when(
      Time < opts_gca$min_time ~ "before",
      opts_gca$max_time < Time ~ "after",
      TRUE ~ "inside")) %>% 
  split(.$window)

ggplot() + 
  aes(x = Time, y = Prop_Target) + 
  hline_chance + 
  vline_onset + 
  stat_mean_se(aes(group = HearsNativeDialect), 
               data = bind_rows(data_windows$before, data_windows$after),
               color = "grey50") + 
  stat_mean_se(aes(color = `Child hears`), data = data_windows$inside) +
  colors +
  ylim(0, 1) +
  theme_ipsum_rc(axis_title_size = 11) + 
  legend_top + 
  labs(x = plot_text$x_time, 
       y = plot_text$y_target) 

ggplot() + 
  aes(x = Time, y = Prop_Target, 
      group = interaction(ResearchID, HearsNativeDialect)) + 
  hline_chance + 
  vline_onset + 
  stat_mean(data = d_m, color = "grey50") + 
  stat_mean(aes(color = `Child hears`), data = data_windows$inside) +
  colors +
  theme_ipsum_rc(axis_title_size = 11) + 
  legend_top + 
  labs(x = plot_text$x_time, 
       y = plot_text$y_target)
```

### Extract orthogonal polynomials

We represent time in our growth curve models as a cubic orthogonal polynomial.
This means that we convert our _Time_ variable into three uncorrelated trend
curves. These are the predictors we use to represent time. I wrote an R package 
for dealing with orthogonal polynomials, and I have a [blog post discussing 
their use in growth curve models](https://tjmahr.github.io/polypoly-package-released/).

```{r, fig.width=4, out.width="50%"}
# Apply analysis window and create orthogonal polynomial of time
d_m <- d_m %>% 
  filter(opts_gca$min < Time, Time < opts_gca$max_time) %>% 
  polypoly::poly_add_columns(Time, 3, prefix = "ot", scale_width = 1)
```

```{r polypoly, fig.width=5, out.width="60%", echo = FALSE}
d_m %>% 
  distinct(Time, ot1, ot2, ot3) %>% 
  tidyr::gather(Degree, value, -Time) %>% 
  ggplot() + 
    aes(x = Time, y = value, color = Degree) + 
    geom_line() + 
    labs(title = "Trial-time features") + 
    theme_ipsum_rc(axis_title_size = 11, plot_title_size = 13) +
    labs(x = plot_text$x_time,
         y = "Feature value")
```

### Model description 

We are going to use a _generalized linear mixed effects model_ to perform the
growth curve analysis. I'll try to summarize what this means very briefly, but
Mirman's 2014 book is the standard reference on the topic.

#### "Generalized linear"

We assume that the data is generated by a binomial process, where the number of
sucesses (looks to target) depends on the total number of attempts (the number
of looks) and the probability of sucess (fixating on the target). This process
assumption shows up in the `glmer()` call as `family = binomial`. The number of
successes and number of attempts show up in the `glmer()` formula with response
variable `cbind(Primary, Others)`: The number of looks to the primary AOI
(target) versus the number of looks to the others.

Our linear model estimates the probability part of the binomial process: At each
time step, what is probability of fixating on the target image? If the linear
trend is positive, then this probability of success increases over the course of
the trial. If the condition effect is negative, then the average probability in
one condition is smaller than in the other. But probabilities have a range from
0 to 1, but linear models run on the range of plus-minus infinity. We need a
transformation to get from probability land to linear land and back again. This
is the log-odds or logit transformation (from probability to linear) and the
logistic transformation (from linear back to probability). These two are the
basis for logistic regression.

```{r logit, fig.width=5, out.width="60%", echo = FALSE}
ggplot(data.frame(x = seq(-10, 10, by = .001))) + 
  aes(x = x) + 
  stat_function(fun = plogis, colour = "red") + 
  theme_ipsum_rc(axis_title_size = 11, plot_title_size = 13) +
  labs(x = "log-odds", y = "probability",
       title = "How log-odds map onto probabilities")
```

Instead of modeling the probability of looking to target, we model the log-odds
of looking to target. This transformation lets us _generalize_ our linear model
into probability land and to model data from a binomial process.

#### "Mixed effects"

I have an [extensive blog
post](https://tjmahr.github.io/plotting-partial-pooling-in-mixed-effects-models/)
illustrating what mixed effects models do, so I'll address their specific role
here. We are going to assume that each child's growth curve is drawn from a
distribution of related growth curves. We are also going to assume that each
child's growth curve in each condition is also drawn from a distribution of
related curves. Each child's growth curve is a combination of:

1. An average growth curve: `... ~ ot1 + ot2 + ot3 + ...`
2. Plus a child's general growth curve features: 
   `... ~ ... + (ot1 + ot2 + ot3 | ResearchID)`
3. Plus a child's specific growth curve features in each condition: 
   `... ~ ... + (ot1 + ot2 + ot3 | HearsOwnDialect:ResearchID)`

The first recognizes that there is an average kind of performance across all
children. The second captures how a child will perform similarly on both
conditions, so it's like a child's _latent_ growth curve shape. The third
captures specific condition level variability, so it's like a child's _observed_
growth curve shape. In the formulas below, the shorthand 
`(ot1 + ot2 + ot3 | ResearchID / HearsNativeDialect)` does both (2) and (3) 
with the x-has-y-nested-in-it operator `x / y`.


## Evaluate the condition effect

First, we need to assess the condition effect. We try four models: 

1. Entirely omit condition.
2. Add condition as a fixed effect, so it shapes the average growth curve.
3. Add condition as a random effect but not as a fixed, so it captures each 
   child's condition-level variability.
4. Add condition as both, so that we capture condition-level variability within 
   children but also estimate the effect of condition on the average growth 
   curve shape.

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

cond_omitted <- glmer(
  cbind(Primary, Others) ~ 
    (ot1 + ot2 + ot3) + 
    (ot1 + ot2 + ot3 | ResearchID),
  data = d_m, 
  control = glmer_controls,
  family = binomial)
summary(cond_omitted)

cond_as_fixed_eff <- glmer(
  cbind(Primary, Others) ~ 
    (ot1 + ot2 + ot3) * HearsNativeDialect + 
    (ot1 + ot2 + ot3 | ResearchID),
  data = d_m,
  control = glmer_controls,
  family = binomial)
summary(cond_as_fixed_eff)

cond_as_random_eff <- m2b <- glmer(
  cbind(Primary, Others) ~ 
    (ot1 + ot2 + ot3) + 
    (ot1 + ot2 + ot3 | ResearchID / HearsNativeDialect),
  data = d_m,
  control = glmer_controls,
  family = binomial)
summary(cond_as_random_eff)

cond_as_both <- glmer(
  cbind(Primary, Others) ~ 
    (ot1 + ot2 + ot3) * HearsNativeDialect + 
    (ot1 + ot2 + ot3 | ResearchID / HearsNativeDialect),
  data = d_m,
  control = glmer_controls,
  family = binomial)
summary(cond_as_both)
```

Use a model comparison on the four variants.

```{r, echo = FALSE}
anova_results <- anova(cond_omitted, cond_as_fixed_eff, 
                       cond_as_random_eff, cond_as_both) %>% 
  as.data.frame() %>% 
  tibble::rownames_to_column("model") %>% 
  as_data_frame()

anova_results %>% 
  select(-deviance) %>% 
  mutate_at(c("AIC", "BIC", "logLik"), round) %>% 
  mutate_at("Chisq", round, 2) %>% 
  mutate_at("Pr(>Chisq)", round, 4) %>% 
  knitr::kable()
```

This is kind of unusual. All of the model comparison metrics favor the model
that includes Child x Native Dialect random effects but omits the Native Dialect
fixed effects. My interpretation:

* For each child, we have two growth curves: The data from hearing their 
  native dialect and the data from hearing the non-native dialect. Their data 
  is _nested_ in these two kinds of blocks.
* This variability is really important for the model, so we include it in 
  the random effects.
* But on average, hearing one's native dialect or not does not influence 
  the growth curve features in a systematic way. So we can ignore using 
  dialect condition as a predictor.
* When we include that native dialect in the model's random effects, we are 
  incorporating information about the nesting of the data. Instead of 
  capturing Child x Native-Dialect effects, it seems like the random effects 
  are capturing Child x Block-of-Testing variability.


## Include child-level predictors

### Data preparation

First, we remove the few participants without EVT or Maternal Education data.

```{r}
to_remove <- d_m %>% 
  filter(is.na(EVT_Standard) | is.na(Maternal_Education_Group)) %>% 
  select(ResearchID, Dialect, Maternal_Education_Group, EVT_Standard) %>% 
  distinct()
to_remove
```

For these models, we are going to have the predictor variables interact with the
intercept term and the linear term, so that the model captures how expressive
vocabulary affects overall accuracy and the slope of the growth curve. We also
scale EVT so that it has a mean of 0 and a standard deviation of 1, so that a
unit change in the EVT predictor represents a 1 standard deviation change in
vocabulary size. We also set maternal education so the mid group is the
reference group.

```{r}
d_narrow <- anti_join(d_m, to_remove) %>% 
  mutate(EVTc = scale(EVT_Standard) %>% as.vector(),
         Medu = factor(Maternal_Education_Group, c("Mid", "Low", "High")))
```

To summarise these measures

```{r}
d_narrow %>% 
  distinct(ResearchID, Medu, EVT_Standard, EVTc, Dialect) %>% 
  group_by(Dialect) %>% 
  summarise(n = n(),
            EVT = mean(EVT_Standard) %>% round(),
            EVT_sd = sd(EVT_Standard) %>% round(),
            EVT_scale = mean(EVTc) %>% round(2),
            EVT_scale_sd = sd(EVTc) %>% round(2),
            Min_EVT = min(EVT_Standard),
            Max_EVT = max(EVT_Standard)) %>% 
  ungroup() %>% 
  knitr::kable()

d_narrow %>% 
  distinct(ResearchID, Medu, EVT_Standard, EVTc, Dialect) %>% 
  group_by(Dialect, Medu) %>% 
  summarise(n = n(),
            EVT = mean(EVT_Standard) %>% round(),
            EVT_sd = sd(EVT_Standard) %>% round(),
            EVT_scale = mean(EVTc) %>% round(2),
            EVT_scale_sd = sd(EVTc) %>% round(2),
            Min_EVT = min(EVT_Standard),
            Max_EVT = max(EVT_Standard)) %>% 
  ungroup() %>% 
  knitr::kable()
```

### EVT and Dialect effects

First we model the effect of EVT.

```{r evt, cache = TRUE}
m_evt <- glmer(
  cbind(Primary, Others) ~ 
    (1 + ot1) * EVTc + ot2 + ot3 +
    (ot1 + ot2 + ot3 | ResearchID / HearsNativeDialect),
  data = d_narrow,
  control = glmer_controls,
  family = binomial)
summary(m_evt)
```

This model tells us that increasing expressive vocabulary by 1 SD improves
overall accuracy and processing efficiency. This is entirely as expected.

```{r evt-preds, echo = FALSE}
new_data_template <- d_narrow %>% 
  distinct(Dialect, Time, ot1, ot2, ot3) %>% 
  mutate(ResearchID = "FakeChild")

predictions <- new_data_template %>% 
  tidyr::crossing(EVTc = -1:1) %>% 
  mutate(fitted = predict(m_evt, newdata = ., re.form = ~ 0, type = "resp"))

ggplot(predictions) + 
  aes(x = Time, y = fitted, linetype = factor(EVTc)) + 
  hline_chance + 
  geom_line(size = 1) + 
  theme_ipsum_rc(axis_title_size = 11, plot_title_size = 13) +
  labs(
    x = plot_text$x_time,
    linetype = "EVT (z score)", 
    y = "Predicted proportion of looks",
    caption = "Conditioned on vocabulary") 
```

Now we include the child's native dialect as a predictor of accuracy and
processing efficieny. We also allow dialect to interact with vocabulary size.
Because the two groups of children differ in the vocabulary levels, we expect
one of the features to be redundant.

```{r dialect, cache = TRUE}
m_dialect <- glmer(
  cbind(Primary, Others) ~ 
    (1 + ot1) * Dialect + ot2 + ot3 +
    (ot1 + ot2 + ot3 | ResearchID / HearsNativeDialect),
  data = d_narrow,
  control = glmer_controls,
  family = binomial)
summary(m_dialect)

m_dialect_evt <- glmer(
  cbind(Primary, Others) ~ 
    (1 + ot1) * Dialect + (1 + ot1) * EVTc + ot2 + ot3 +
    (ot1 + ot2 + ot3 | ResearchID / HearsNativeDialect),
  data = d_narrow,
  control = glmer_controls,
  family = binomial)
summary(m_dialect_evt)

m_dialect_x_evt <- glmer(
  cbind(Primary, Others) ~ 
    (1 + ot1) * Dialect * EVTc + ot2 + ot3 +
    (ot1 + ot2 + ot3 | ResearchID / HearsNativeDialect),
  data = d_narrow,
  control = glmer_controls,
  family = binomial)
summary(m_dialect_x_evt)
```

It's hard to tell what is going on without plotting. First, look at the dialect 
only effects.

```{r dialect-effects, echo = FALSE}
predictions <- new_data_template %>% 
  tidyr::crossing(EVTc = -1:1) %>% 
  mutate(fitted = predict(m_dialect, newdata = ., re.form = ~ 0, type = "resp"))

ggplot(predictions) + 
  aes(x = Time, y = fitted, color = Dialect) + 
  hline_chance + 
  geom_line(size = 1) + 
  theme_ipsum_rc(axis_title_size = 11, plot_title_size = 13) +
  labs(
    x = plot_text$x_time,
    y = "Predicted proportion of looks",
    caption = "Conditioned on native dialect") 

predictions <- new_data_template %>% 
  tidyr::crossing(EVTc = -1:1) %>% 
  mutate(fitted = predict(m_dialect_evt, newdata = ., 
                          re.form = ~ 0, type = "resp"))

ggplot(predictions) + 
  aes(x = Time, y = fitted, color = Dialect, linetype = factor(EVTc)) + 
  hline_chance + 
  geom_line(size = 1) + 
  theme_ipsum_rc(axis_title_size = 11, plot_title_size = 13) +
  labs(
    x = plot_text$x_time,
    linetype = "EVT (z score)", 
    y = "Predicted proportion of looks",
    caption = "Conditioned on native dialect and vocabulary") 

predictions <- new_data_template %>% 
  tidyr::crossing(EVTc = -1:1) %>% 
  mutate(fitted = predict(m_dialect_x_evt, newdata = ., 
                          re.form = ~ 0, type = "resp"))

ggplot(predictions) + 
  aes(x = Time, y = fitted, color = Dialect, linetype = factor(EVTc)) + 
  hline_chance + 
  geom_line(size = 1) + 
  theme_ipsum_rc(axis_title_size = 11, plot_title_size = 13) +
  labs(
    x = plot_text$x_time,
    linetype = "EVT (z score)", 
    y = "Predicted proportion of looks",
    caption = "Conditioned on native dialect x vocabulary") 
```

This final plot would indicate that the MAE group performs very similarly, so
that vocabulary size does not adequately differentiate these children, whereas a 1
SD change in the AAE group predicts a large improvement in word recognition.
However, we need to take these predictions with a grain of salt because the -1
z-EVT MAE speakers and the +1 z-EVT AAE speakers are not represented in the
data.


#### Model comparison 

Model comparison indicates the dialect is a better child level predictor than 
expressive vocabulary, and that including simple effects of dialect and
expressive vocabulary is redundant (which we expected).

```{r, echo = FALSE}
tidy_anova <- . %>% 
  as.data.frame() %>% 
  tibble::rownames_to_column("model") %>% 
  as_data_frame() %>% 
  select(-deviance) %>% 
  mutate_at(c("AIC", "BIC", "logLik"), round) %>% 
  mutate_at("Chisq", round, 2) %>% 
  mutate_at("Pr(>Chisq)", round, 4) %>% 
  knitr::kable()
  
anova(m_evt, m_dialect, m_dialect_evt, m_dialect_x_evt) %>% 
  tidy_anova()
```

Two of the three model comparison metrics, however, prefer the model with the
interaction effect.

```{r, echo = FALSE}
anova(m_evt, m_dialect, m_dialect_x_evt) %>% 
  tidy_anova()
```

### Include Maternal education

Now we include maternal education as a predictor and compare it against dialect.

```{r medu, cache = TRUE}
m_medu <- glmer(
  cbind(Primary, Others) ~ 
    (1 + ot1) * Medu + ot2 + ot3 +
    (ot1 + ot2 + ot3 | ResearchID / HearsNativeDialect),
  data = d_narrow,
  control = glmer_controls,
  family = binomial)
summary(m_medu)

m_medu_dialect <- glmer(
  cbind(Primary, Others) ~ 
    (1 + ot1) * Medu + (1 + ot1) * Dialect + 
    ot2 + ot3 +
    (ot1 + ot2 + ot3 | ResearchID / HearsNativeDialect),
  data = d_narrow,
  control = glmer_controls,
  family = binomial)
summary(m_medu_dialect)

m_medu_x_dialect <- glmer(
  cbind(Primary, Others) ~ 
    (1 + ot1) * Medu * Dialect + 
    ot2 + ot3 +
    (ot1 + ot2 + ot3 | ResearchID / HearsNativeDialect),
  data = d_narrow,
  control = glmer_controls,
  family = binomial)
summary(m_medu_x_dialect)
```


#### Model comparison 

Model comparison favors the model _without_ the dialect by maternal education
effect.

```{r, echo = FALSE}
anova(m_dialect, m_medu, m_medu_dialect, m_medu_x_dialect) %>% tidy_anova()
anova(m_dialect, m_medu, m_medu_x_dialect) %>% tidy_anova()
```

#### Model fits

As before, the models are best understood through plotting.

```{r medu-preds, echo = FALSE}
new_data_template <- d_narrow %>% 
  distinct(Dialect, Medu, Time, ot1, ot2, ot3) %>% 
  mutate(ResearchID = "FakeChild") %>% 
  mutate(MeduPlot = factor(Medu, c("Low", "Mid", "High")))

predictions <- new_data_template %>% 
  tidyr::crossing(EVTc = -1:1) %>% 
  mutate(fitted = predict(m_medu, newdata = ., re.form = ~ 0, type = "resp"))

ggplot(predictions) + 
  aes(x = Time, y = fitted, color = MeduPlot) + 
  hline_chance + 
  geom_line(size = 1) + 
  theme_ipsum_rc(axis_title_size = 11, plot_title_size = 13) +
  viridis::scale_color_viridis(discrete = TRUE, end = .7, option = "C") +
  labs(
    x = plot_text$x_time,
    color = "Maternal Ed.", 
    y = "Predicted proportion of looks",
    caption = "Conditioned on maternal education") 

predictions <- new_data_template %>% 
  tidyr::crossing(EVTc = -1:1) %>% 
  mutate(fitted = predict(m_medu_dialect, newdata = ., re.form = ~ 0, type = "resp"))

ggplot(predictions) + 
  aes(x = Time, y = fitted, color = MeduPlot, linetype = Dialect) + 
  hline_chance + 
  geom_line(size = 1) + 
  theme_ipsum_rc(axis_title_size = 11, plot_title_size = 13) +
  viridis::scale_color_viridis(discrete = TRUE, end = .7, option = "C") +
  labs(
    x = plot_text$x_time,
    color = "Maternal Ed.", 
    y = "Predicted proportion of looks",
    caption = "Conditioned on maternal ed. and native dialect") 
```

Some implications from these models:

* There is no difference between mid and high maternal education. 
* The dialect effect disappears when we account for maternal education.

## Maternal education

Finally, we pit expressive vocabulary against maternal education. 

```{r evt-dialect, cache = TRUE}
m_medu_evt <- glmer(
  cbind(Primary, Others) ~ 
    (1 + ot1) * Medu + (1 + ot1) * EVTc + 
    ot2 + ot3 +
    (ot1 + ot2 + ot3 | ResearchID / HearsNativeDialect),
  data = d_narrow,
  control = glmer_controls,
  family = binomial)
summary(m_medu_evt)

m_medu_x_evt <- glmer(
  cbind(Primary, Others) ~ 
    (1 + ot1) * Medu * EVTc + 
    ot2 + ot3 +
    (ot1 + ot2 + ot3 | ResearchID / HearsNativeDialect),
  data = d_narrow,
  control = glmer_controls,
  family = binomial)
summary(m_medu_x_evt)
```

None of the effects of EVT appear to model from the model summary. Model
comparison provides no support for the model with EVT and maternal education.
Only one of the three model comparison metrics support the model with EVT and
maternal education.

```{r, echo = FALSE}
anova(m_medu, m_medu_evt, m_medu_x_evt) %>% tidy_anova()
anova(m_dialect, m_medu_x_evt) %>% tidy_anova()
```

It might be worth "refitting" the last model with low group as the reference
condition to see what happens. (The results or predictions don't change... just
the way the values are reported changed.) This model shows that the expressive
vocabulary x Time effect is significant in the low maternal education group.

```{r medu2, cache = TRUE}
d_narrow2 <- d_narrow %>% 
  mutate(Medu = factor(Medu, c("Low", "Mid", "High")))

m_medu_x_evt2 <- glmer(
  cbind(Primary, Others) ~ 
    (1 + ot1) * Medu * EVTc + 
    ot2 + ot3 +
    (ot1 + ot2 + ot3 | ResearchID / HearsNativeDialect),
  data = d_narrow2,
  control = glmer_controls,
  family = binomial)
summary(m_medu_x_evt2)
```


Again, we plot the model fits.

```{r}
predictions <- new_data_template %>% 
  tidyr::crossing(EVTc = -1:1) %>% 
  mutate(fitted = predict(m_medu_evt, newdata = ., re.form = ~ 0, type = "resp"))

ggplot(predictions) + 
  aes(x = Time, y = fitted, color = MeduPlot, linetype = factor(EVTc)) + 
  hline_chance + 
  geom_line(size = 1) + 
  viridis::scale_color_viridis(discrete = TRUE, end = .7, option = "C") +
  theme_ipsum_rc(axis_title_size = 11, plot_title_size = 13) +
  labs(
    x = plot_text$x_time,
    linetype = "EVT (z-score)",
    color = "Maternal Ed.", 
    y = "Predicted proportion of looks",
    caption = "Conditioned on maternal ed. and vocabulary") 

predictions <- new_data_template %>% 
  tidyr::crossing(EVTc = -1:1) %>% 
  mutate(fitted = predict(m_medu_x_evt, newdata = ., 
                          re.form = ~ 0, type = "resp"))
ggplot(predictions) + 
  aes(x = Time, y = fitted, color = MeduPlot, linetype = factor(EVTc)) + 
  hline_chance + 
  geom_line(size = 1) + 
  viridis::scale_color_viridis(discrete = TRUE, end = .7, option = "C") +
  theme_ipsum_rc(axis_title_size = 11, plot_title_size = 13) +
  labs(
    x = plot_text$x_time,
    linetype = "EVT (z-score)",
    color = "Maternal Ed.", 
    y = "Predicted proportion of looks",
    caption = "Conditioned on maternal ed. x vocabulary") 
```

Here's what the interaction between maternal education and expressive vocabulary
is capturing: The blue lines (the high maternal education group) are all
complete on top of each other. There appears to be no differentiating effects of
vocabulary in this group. But for the low maternal education group, expressive
vocabulary does make a difference.

## Findings

* No effect of experiment condition.
* Our dialect groups differed in maternal education and vocabulary size too.
  Each of these three measures predicted growth curve features but the best 
  single predictor among them is maternal education.
* Vocabulary size might matter more the low maternal education group but the 
  evidence is not overwhelming.


***

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