---
title: "Observing interviewer performance in slices or by traces: A comparison of methods to predict interviewers’ individual contributions to interviewer variance"
date: "`r format(Sys.time(), '%d %B %Y')`"
author: "Celine Wuyts and Geert Loosveldt"
geometry: margin = 3cm
output:
  rmdformats::downcute:
    number_sections: TRUE
    df_print: paged
---

```{r setup, include = FALSE}

knitr::opts_chunk$set(error = FALSE, warning = FALSE, message = FALSE,
                      fig.width = 7, fig.height = 10)

library(tidyverse)
library(ggthemes)
library(cowplot)
library(stargazer)
library(formattable)
library(kableExtra)

library(lme4)
library(merTools)

library(sjstats)
library(broom)

options(knitr.kable.NA = '.')

```

```{r themes, include = FALSE}

themeKUL <- theme_tufte(base_size = 9) + # base_family = "Calibri"
    theme(axis.title = element_text(size = 9, face = "plain"),
          axis.text = element_text(size = 9),
          axis.line.x = element_line(),
          plot.title = element_blank(),
          legend.title = element_blank(),
          legend.text = element_text(size = 9),
          strip.text = element_text(size = 9, face = "bold"),
          legend.position = "bottom",
          legend.direction = "horizontal",
          legend.box = "vertical",
          legend.spacing = unit(0, "line"),
          legend.key.size = unit(.75, "line"))

```


# Data and methods

```{r get_data}

# ICCs
ICCs <- read.csv("20200108_ICCs.csv", dec = ".", 
                 sep = ";", stringsAsFactors = FALSE)

ICCs <- ICCs %>%
  filter(nobs >= 300) 
# Drops 13 items that are routed to small groups of resondents

# Random effects
RE <- read.csv("20200108_RE.csv", dec = ".", 
                 sep = ";", stringsAsFactors = FALSE)

# Interviewers
interviewers <- read.csv("20200226_interviewers.csv" , dec = ".", 
                 sep = ";", stringsAsFactors = FALSE)


# Cumulative interview time measures
speed <- read.csv("speedindicators.ICum.csv" , dec = ".", 
                 sep = ";", stringsAsFactors = FALSE)

RE_cum <- merge(merge(RE, interviewers[c("essround", "cntry", "intnum", "INTNUM_ANON", "iWL")], 
                      by = c("essround", "cntry", "intnum")),
                speed, by = c("essround", "cntry", "intnum"), all.x = T)
```

```{r audio_alternative groupings}

interviewers <- interviewers %>%
  mutate(AUDIO_READING = ((AUDIO_ALL == 2) +
           (AUDIO_ALLQSAME == 2) +
           (AUDIO_ALLREAD == 2) +
           (AUDIO_CLEAR == 2) +
           (AUDIO_COMPLETE == 2) + 
           (AUDIO_NOEXTRA == 2)),
         AUDIO_SUGGESTIVE = ((AUDIO_NOEXAMPLE == 2) +
                               (AUDIO_NOJUDGEMENT == 2) +
                               (AUDIO_NOOPINION == 2) +
                               (AUDIO_NOTINOPTIONS == 2) +
                               (AUDIO_QREPEAT == 2) +
                               (AUDIO_NOTSUGGESTIVE == 2) +
                               (AUDIO_RESPINTERPRET == 2)),
         AUDIO_RUSHING = ((AUDIO_TEMPO == 2) +
                            (AUDIO_ENOUGHTIME == 2)),
         AUDIO_READINGCAT = ifelse(AUDIO_READING > 0, 1, 0),
         AUDIO_SUGGESTIVECAT = ifelse(AUDIO_SUGGESTIVE > 0, 1, 0),
         AUDIO_RUSHINGCAT = ifelse(AUDIO_RUSHING > 0, 1, 0))

```

## Selection of survey items

```{r select_items}

# Spread by round
ICCs_itembyround <- ICCs %>%
  mutate(essround = paste0("ESS", essround, "_BEDUT")) %>%
  dplyr::select(essround, item, icc_fit_controls, label) %>%
  spread(key = essround, value = icc_fit_controls)

# Number of items
with(ICCs_itembyround, addmargins(table(!is.na(ESS6_BEDUT), !is.na(ESS7_BEDUT))))

# Subset items ICC > 3%
ICCs_selected <- ICCs %>%
  filter(LRtest_pvalue < 0.05 & icc_fit_controls > 0.03)   

# Spread by round
ICCs_selected_itembyround <- ICCs_selected %>%
  mutate(essround = paste0("ESS", essround, "_BEDUT")) %>%
  dplyr::select(essround, item, icc_fit_controls, label) %>%
  spread(key = essround, value = icc_fit_controls)

# Number of items selected
with(ICCs_selected_itembyround, addmargins(table(!is.na(ESS6_BEDUT), !is.na(ESS7_BEDUT))))

```

```{r ICCs_overviewtable, results = "asis"}

overviewtable <- ICCs_selected_itembyround %>%
  mutate(icc_fit_controls_mean = rowMeans(cbind(ESS6_BEDUT, ESS7_BEDUT), na.rm = T)) %>% 
  arrange(-icc_fit_controls_mean) %>%
  dplyr::select(item, label, ESS6_BEDUT, ESS7_BEDUT) %>%
  mutate(label = gsub("/", ", ", label))

kable(overviewtable,
      row.names = FALSE,
      escape = TRUE,
      digits = 3,
      col.names = c("Name", 
                    "Label",
                    "ESS6 BEDUT",
                    "ESS7 BEDUT"),
      align = c(rep("l", 2), rep("r", 2)),
      caption = "Intra-interviewer correlations for selected survey items") %>%
  column_spec(c(1, 3:4), width = "1.5cm") %>%
  column_spec(2, width = "8cm") %>%
  add_header_above(c(" " = 2, "Intra-interviewer correlation" = 2))

```

```{r ICCs_descriptives, results = "asis"}

stargazer(ICCs["icc_fit_controls"],
          type = "html",
          header = FALSE,
          median = TRUE, iqr = TRUE,
          title = "Descriptives of intra-interviewer correlations (all survey items)")
```

```{r ICCs_descriptives_byround, results = "asis"}

stargazer(ICCs_itembyround[c("ESS6_BEDUT", "ESS7_BEDUT")],
          type = "html",
          header = FALSE,
          median = TRUE, iqr = TRUE,
          title = "Descriptives of intra-interviewer correlations (all survey items)")
```

```{r ICCs_selected_descriptives, results = "asis"}

stargazer(ICCs_selected["icc_fit_controls"],
          type = "html",
          header = FALSE,
          median = TRUE, iqr = TRUE,
          title = "Descriptives of intra-interviewer correlations (selected survey items)")
```

```{r ICCs_selected_descriptives_byround, results = "asis"}

stargazer(ICCs_selected_itembyround[c("ESS6_BEDUT", "ESS7_BEDUT")],
          type = "html",
          header = FALSE,
          median = TRUE, iqr = TRUE,
          title = "Descriptives of intra-interviewer correlations (selected survey items)")
```


```{r structure_data_interviewers}

audioscore_byround <- interviewers %>%
  mutate(essround = paste0("ESS", essround, "_BEDUT")) %>%
  dplyr::select(essround, INTNUM_ANON, AUDIO_NDEV) %>%
  spread(key = essround, value = AUDIO_NDEV)

speed_byround <- interviewers %>%
  mutate(essround = paste0("ESS", essround, "_BEDUT")) %>%
  dplyr::select(essround, INTNUM_ANON, speed_mainexclF.mean) %>%
  spread(key = essround, value = speed_mainexclF.mean)

speedtrend_byround <- interviewers %>%
  mutate(essround = paste0("ESS", essround, "_BEDUT")) %>%
  dplyr::select(essround, INTNUM_ANON, speedTrend.mean) %>%
  spread(key = essround, value = speedTrend.mean)

```

## Interviewers' deviation from the standardized interviewing protocol

```{r audio_checks_frequencies}

audiochecks_freq <- read.csv("20200108_audiochecks_freq.csv", dec = ".", 
                 sep = ";", stringsAsFactors = FALSE)

audiovarlabels <- read.csv("var_audio labels.csv", dec = ".", 
                 sep = ";", stringsAsFactors = FALSE)

audiovarlabels <- audiovarlabels%>%
  mutate(checkorder = 1:nrow(audiovarlabels))

audiochecks_freq_checksbyround <- audiochecks_freq %>%
  mutate(essround = paste0("ESS", essround, "_BEDUT")) %>%
  dplyr::select(essround, variable, label, prop) %>%
  spread(key = essround, value = prop)

audiochecks_freq_checksbyround <- left_join(audiochecks_freq_checksbyround, audiovarlabels) %>%
  arrange(checkorder)

```

```{r audio_freq, results = "asis"}

audiotable <- audiochecks_freq_checksbyround %>%
        dplyr::select(label, ESS6_BEDUT, ESS7_BEDUT) 

kable(audiotable %>%
        mutate_at(c("ESS6_BEDUT", "ESS7_BEDUT"), 
                  function(x) formattable::percent(x, digits = 1)),
      row.names = FALSE,
      booktabs = TRUE,
      digits = 1,
      col.names = c("Check",
                    "ESS6 BEDUT",
                    "ESS7 BEDUT"),
      align = c("l", rep("r", 2)),
      caption = "Checklist of standardized interviewing protocol") %>%
  column_spec(1, width = "10cm") %>%
  column_spec(2:3, width = "2.5cm") %>%
  add_header_above(c(" " = 1, "Deviation from check" = 2), escape = FALSE)

```

```{r audio_descriptives, results = "asis"}

stargazer(audioscore_byround[c("ESS6_BEDUT", "ESS7_BEDUT")],
          type = "html",
          header = FALSE,
          title = "Descriptives of standardized interviewing deviation score")
```

```{r audio_counts}

with(interviewers, prop.table(table(AUDIO_NDEV, essround), 2))

with(interviewers, prop.table(table(AUDIO_READING > 0, essround), 2))
with(interviewers, prop.table(table(AUDIO_SUGGESTIVE > 0, essround), 2))
with(interviewers, prop.table(table(AUDIO_RUSHING > 0, essround), 2))

```

```{r audio_histogram, fig.cap = "Distribution of the interviewers' standardized interviewing deviation score", fig.height = 3}

p <- ggplot(data = interviewers %>%
              mutate(essround = as.factor(paste0("ESS", essround, "_BEDUT"))),
            aes(x = AUDIO_NDEV)) +
  geom_bar(alpha = .5, col = NA, fill = "grey", width = 0.75) +
  geom_hline(yintercept = c(5, 10, 15), colour = 'white', width = 0.25) +
  scale_x_continuous(name = "Standardized interviewing deviation score", 
                     breaks = seq(0, 30, 1)) +
  facet_wrap(facets = "essround") +
  themeKUL +
  theme(axis.ticks.y = element_blank(),
        axis.title.y = element_blank(),
        axis.line.x = element_line(),
        legend.position = "none")
p

```

```{r audio_scatterplot, fig.cap = "Joint distribution of the two survey rounds' standardized interviewing deviation score", fig.height = 4, fig.width = 4}

p <- ggplot(data = subset(audioscore_byround, complete.cases(audioscore_byround)),
              aes(x = ESS6_BEDUT, y = ESS7_BEDUT)) +
  geom_point(size = 2, shape = 21, alpha = .2, fill = "black", col = "white") + 
  scale_x_continuous(name = "Standardized interviewing deviation score, ESS6_BEDUT", 
                     limits = c(-1, 15), breaks = seq(0, 15, 1)) +
  scale_y_continuous(name = "Standardized interviewing deviation score, ESS7_BEDUT", 
                     limits = c(-1, 15), breaks = seq(0, 15, 1)) +
  themeKUL +
  theme(axis.line.x = element_line(),
        axis.line.y = element_line()) 

xplot <- ggplot(subset(audioscore_byround, complete.cases(audioscore_byround)),
                aes(x = ESS6_BEDUT)) + 
  geom_histogram(col = "white", alpha = .5, binwidth = 1) + 
  scale_x_continuous(limits = c(-1, 15)) + 
  themeKUL + 
  theme(axis.line = element_blank(),
        axis.ticks = element_blank(),
        axis.title = element_blank(),
        axis.text = element_blank())

yplot <- ggplot(subset(audioscore_byround, complete.cases(audioscore_byround)),
                aes(x = ESS7_BEDUT)) + 
  geom_histogram(col = "white", alpha = .5, binwidth = 1) +
  scale_x_continuous(limits = c(-1, 15)) + 
  themeKUL + 
  theme(axis.line = element_blank(),
        axis.ticks = element_blank(),
        axis.title = element_blank(),
        axis.text = element_blank()) + 
  coord_flip()

p <- plot_grid(xplot, NULL, p, yplot, ncol = 2, align = "hv", 
      rel_widths = c(3, 1), rel_heights = c(1, 3))

p

```


## Interviewers' interview speed

```{r speed_descriptives, results = "asis"}

stargazer(speed_byround[c("ESS6_BEDUT", "ESS7_BEDUT")],
          type = "html",
          header = FALSE,
          title = "Descriptives of average interview speed")
```

```{r speedtrend_descriptives, results = "asis"}

stargazer(speedtrend_byround[c("ESS6_BEDUT", "ESS7_BEDUT")],
          type = "html",
          header = FALSE,
          title = "Descriptives of average interview acceleration over the questionnaire")
```

```{r speed_scatterplot, fig.cap = "Joint distribution of the two survey rounds' average interview speed", fig.height = 4, fig.width = 4}

p <- ggplot(data = subset(speed_byround, complete.cases(speed_byround)),
              aes(x = ESS6_BEDUT, y = ESS7_BEDUT)) +
  geom_point(size = 2, shape = 21, alpha = .2, fill = "black", col = "white") + 
  scale_x_continuous(name = "Average interview speed, ESS6_BEDUT", 
                     limits = c(0, 6), breaks = seq(0, 6, 1)) +
  scale_y_continuous(name = "Average interview speed, ESS7_BEDUT", 
                     limits = c(0, 6), breaks = seq(0, 6, 1)) +
  themeKUL +
  theme(axis.line.x = element_line(),
        axis.line.y = element_line()) 

xplot <- ggplot(subset(speed_byround, complete.cases(speed_byround)),
                aes(x = ESS6_BEDUT)) + 
  geom_histogram(col = "white", alpha = .5, binwidth = .5) + 
  scale_x_continuous(limits = c(0, 6)) + 
  themeKUL + 
  theme(axis.line = element_blank(),
        axis.ticks = element_blank(),
        axis.title = element_blank(),
        axis.text = element_blank())

yplot <- ggplot(subset(speed_byround, complete.cases(speed_byround)),
                aes(x = ESS7_BEDUT)) + 
  geom_histogram(col = "white", alpha = .5, binwidth = .5) +
  scale_x_continuous(limits = c(0, 6)) + 
  themeKUL + 
  theme(axis.line = element_blank(),
        axis.ticks = element_blank(),
        axis.title = element_blank(),
        axis.text = element_blank()) + 
  coord_flip()

p <- plot_grid(xplot, NULL, p, yplot, ncol = 2, align = "hv", 
      rel_widths = c(3, 1), rel_heights = c(1, 3))

p

```

```{r stabilization_speed}

p <- ggplot(data = subset(speed, intnum %in% c(111166, 111247, 112734)) %>%
               mutate(intnum = as.factor(intnum)),
             aes(x = intervieworder, y = speed_mainexclF.mean, linetype = intnum)) +
  geom_line() +
  scale_x_continuous(name = "Number of completed interviews", 
                     limits = c(0, 20), 
                     breaks = 0:20,
                     labels = c(0:10, rep("", 4), 15, rep("", 4), 20)) +
  scale_y_continuous(name = "Average interview speed", limits = c(2, NA), breaks = seq(0, 5, .5)) +
  scale_linetype_manual(values = c("111166" = "dashed",
                                   "111247" = "solid",
                                   "112734" = "dotted"),
                     labels = c("111166" = "Interviewer 1",
                                "111247" = "Interviewer 2",
                                "112734" = "Interviewer 3")) +
  themeKUL +
  theme(axis.line = element_line(),
        legend.direction = "vertical")

p

```

```{r stabilization_acceleration}

p <- ggplot(data = subset(speed, intnum %in% c(111166, 111247, 112734)) %>%
               mutate(intnum = as.factor(intnum)),
             aes(x = intervieworder, y = speedTrend.mean, linetype = intnum)) +
  geom_line() +
  scale_x_continuous(name = "Number of completed interviews", 
                     limits = c(0, 20), 
                     breaks = 0:20,
                     labels = c(0:10, rep("", 4), 15, rep("", 4), 20)) +
  scale_y_continuous(name = "Average acceleration over questionnaire", breaks = seq(-.5, .5, .25)) +
  scale_linetype_manual(values = c("111166" = "dashed",
                                   "111247" = "solid",
                                   "112734" = "dotted"),
                     labels = c("111166" = "Interviewer 1",
                                "111247" = "Interviewer 2",
                                "112734" = "Interviewer 3")) +
  themeKUL +
  theme(axis.line = element_line(),
        legend.direction = "vertical")

p

```


## Interviewers' contributions to measurement error

```{r attach_ICCs}

RE <- left_join(RE, interviewers, 
                by = c("essround", "cntry", "intnum"))
RE <- left_join(RE, ICCs[c("essround", "cntry", "item", "icc_fit_controls")],
                 by = c("essround", "cntry", "item"))

```

```{r trimmed_values}

# Deviations larger than 0.5 in absolute value

with(RE, table(abs_condval > .5))
with(RE, prop.table(table(abs_condval > .5)))

```

```{r RE_plot, fig.cap = "Interviewers' individual absolute bias estimates"}

p <- ggplot(RE %>%
               mutate(abs_condval_mean = ave(abs_condval, INTNUM_ANON, FUN = mean)) %>%
               arrange(abs_condval_mean) %>%
               mutate(essround = as.factor(paste0("ESS", essround, "_BEDUT")),
                      INTNUM_ANON = factor(as.factor(INTNUM_ANON), levels = unique(INTNUM_ANON)),
                      abs_condval = ifelse(abs_condval > .5, .5, abs_condval)), 
            aes(x = abs_condval, y = INTNUM_ANON, col = essround, shape = essround)) +
  geom_point(alpha = .25, col = "black") +
  scale_x_continuous(name = "Estimated absolute difference between interviewer-specific intercept and
                     overall intercept", limits = c(0, .51), breaks = seq(0, .5, .05)) +
  scale_shape_manual(values = c("ESS6_BEDUT" = 2, "ESS7_BEDUT" = 3)) +
  geom_point(inherit.aes = FALSE, aes(x = abs_condval_mean, y = INTNUM_ANON), shape = 8, size = 1, col = "black") +
  themeKUL +
  theme(axis.ticks.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        axis.line.x = element_line(),
        strip.text = element_blank(),
        legend.position = "bottom",
        legend.direction = "vertical") +
  guides(colour = guide_legend(override.aes = list(alpha = .5)),
         shape = guide_legend(override.aes = list(alpha = .5)))

p

```
```{r standardization}

RE_scaled <- RE %>%
  mutate_if(is.numeric, scale)

```

# Results

```{r models}

# Model 0
spec_0 <- "abs_condval ~ essround + iWL + (1|item) + (1|INTNUM_ANON)"

fit_0 <- lmer(spec_0, data = RE_scaled, na.action = na.exclude)
fit_00 <- update(fit_0, . ~ . - iWL)

# Model A and A'
fit_audio <- update(fit_0, . ~ . + AUDIO_NDEV)
fit_audioalt <- update(fit_0, . ~ . + AUDIO_READINGCAT + AUDIO_SUGGESTIVECAT + AUDIO_RUSHINGCAT)

# Model B
fit_speed <- update(fit_0, . ~ . + speed_mainexclF.mean + speedTrend.mean)

anova(fit_audio, fit_0)
anova(fit_speed, fit_0)
anova(fit_audio, fit_speed)

# Model C and C'
fit_final <- update(fit_0, . ~ . + AUDIO_NDEV + speed_mainexclF.mean + speedTrend.mean)
fit_finalalt <- update(fit_0, . ~ . + AUDIO_READINGCAT + AUDIO_SUGGESTIVECAT + AUDIO_RUSHINGCAT + speed_mainexclF.mean + speedTrend.mean)

anova(fit_final, fit_audio)
anova(fit_final, fit_speed)

# # Additional interaction term tests
# fit_speed_quadratic <- update(fit_audio, . ~ . + I(speed_mainexclF.mean^2) + I(speedTrend.mean^2))
# fit_audio_interact <- update(fit_0, . ~ . + AUDIO_NDEV*iWL + AUDIO_NDEV*icc_fit_controls)
# fit_speed_interact <- update(fit_0, . ~ . + speed_mainexclF.mean*iWL + speedTrend.mean*iWL +
#                       speed_mainexclF.mean*icc_fit_controls + speedTrend.mean*icc_fit_controls)

# # Add workload interaction
# fit_audio3 <- update(fit_0, . ~ . + iWL*AUDIO_NDEV)
# fit_speed3 <- update(fit_0, . ~ . + iWL*speed_mainexclF.mean + iWL*speedTrend.mean)
# fit_speedtrend3 <- update(fit_0, . ~ . + iWL*speedTrend.mean)
# fit_final3 <- update(fit_0, . ~ . + iWL*AUDIO_NDEV + iWL*speed_mainexclF.mean + iWL*speedTrend.mean)
# fit_finalalt3 <- update(fit_0, . ~ . + iWL*AUDIO_READINGCAT + iWL*AUDIO_SUGGESTIVECAT + iWL*AUDIO_RUSHINGCAT + iWL*speed_mainexclF.mean + iWL*speedTrend.mean)

# List of model fits
fits <- list("fit_00" = fit_00, "fit_0" = fit_0, 
             "fit_audio" = fit_audio, "fit_audioalt" = fit_audioalt, 
             "fit_speed" = fit_speed, 
             "fit_final" = fit_final, "fit_finalalt" = fit_finalalt)


fits_estimates <- purrr::map_df(fits, tidy, .id = "model")
fits_modelfit <- purrr::map_df(fits, function(f) broom::glance(refitML(f)), .id = "model")
fits_nobs <- purrr::map_df(fits, function(f){
  data.frame(n = nobs(f), J = ngrps(f)[1], K = ngrps(f)[2])}, .id = "model")

# Extract variance components
fits_varcomp <- lapply(fits, function(fit) as.data.frame(VarCorr(fit)))
fits_var_interviewers <- as.character(digits(lapply(fits_varcomp, function(x) x$vcov[x$grp == "INTNUM_ANON"]), 6))
fits_var_items <- as.character(digits(lapply(fits_varcomp, function(x) x$vcov[x$grp == "item"]), 6))
fits_var_residual <- as.character(digits(lapply(fits_varcomp, function(x) x$vcov[x$grp == "Residual"]), 6))
fits_icc_interviewers <- as.character(digits(lapply(fits_varcomp, function(x) x$vcov[x$grp == "INTNUM_ANON"]/sum(x$vcov)), 6))
fits_icc_items <- as.character(digits(lapply(fits_varcomp, function(x) x$vcov[x$grp == "item"]/sum(x$vcov)), 6))

intvar_0 <- fits_varcomp[["fit_0"]]$vcov[fits_varcomp[["fit_0"]]$grp == "INTNUM_ANON"]
fits_intvar_explained <- lapply(fits_varcomp, function(x) (intvar_0 - x$vcov[x$grp == "INTNUM_ANON"])/intvar_0)

# Model fits
fits_AIC <- as.character(comma(lapply(fits, function(fit) AIC(refitML(fit))), digits = 0))

```

```{r models_table, results = "asis"}

stargazer(fit_0, fit_audio, fit_audioalt, fit_speed, fit_final, fit_finalalt,
          type = "html",
          header = FALSE,
          digits = 4,
          title = "Fixed parameter estimates for models explaining interviewers' contributions to interviewer variance",
          model.numbers = FALSE,
          column.labels = c("Model 0", "Model A", "Model A'", "Model B", "Model C", "Model C'"),
          dep.var.caption = "",
          dep.var.labels.include = FALSE,
          order = c("AUDIO_NDEV",
                    "AUDIO_READINGCAT", "AUDIO_SUGGESTIVECAT", "AUDIO_RUSHINGCAT",
                    "speed_mainexclF.mean", "speedTrend.mean",
                    "iWL", "essround", "(Intercept)"),
          covariate.labels = c("AUDIO_NDEV" = "Standardized interviewing deviation count",
                               "AUDIO_READINGCAT" = "At least one flag for question reading",
                               "AUDIO_SUGGESTIVECAT" = "At least one flag for deviation from neutrality",
                               "AUDIO_RUSHINGCAT" = "At least one flag for inadequate interview pace",
                               "speed_mainexclF.mean" = "Average interview speed",
                               "speedTrend.mean" = "Average acceleration over questionnaire",
                               "iWL" = "Number of interviews completed by interviewer",
                               "essround7" = "Round 7"),
          keep.stat = "n",
          add.lines = list(c("Interviewer variance", tail(unlist(fits_var_interviewers), -1)),
                           c("Item variance", tail(unlist(fits_var_items), -1)),
                           c("Residual variance", tail(unlist(fits_var_residual), -1)),
                           c("AIC", tail(unlist(fits_AIC), -1), 0)))
```

```{r cumulativemodels}

AICs_speedmodels <- data.frame()

for(i in 1:max(speed$intervieworder)){
  fit <- lme4::lmer(abs_condval ~ essround + pmin(i, iWL) + (1|item) + (1|INTNUM_ANON) + 
                      speed_mainexclF.mean + speedTrend.mean,
              data = subset(RE_cum, intervieworder == i & iWL >= 8),
              na.action = na.exclude)
    
  assign(paste0("fit_speed_", i), fit)
  AICs_speedmodels <- rbind(AICs_speedmodels, data.frame(intervieworder = i, 
                                                         AIC = AIC(refitML(fit))))
}

fit_0_8 <- lmer(abs_condval ~ essround + iWL + (1|item) + (1|INTNUM_ANON),
              data = subset(RE, iWL >= 8),
              na.action = na.exclude)

AICs_speedmodels_diff <- AICs_speedmodels %>%
  mutate(AICdiff0 = AIC - AIC(refitML(fit_0_8)),
         AICdiffend = AIC - with(subset(AICs_speedmodels, intervieworder == 47), AIC))

```

```{r cumulativemodels_plot, fig.width = 8, fig.height = 6}

p <- ggplot(data = subset(AICs_speedmodels_diff, intervieworder <= 8), 
            aes(y = AICdiff0, x = intervieworder)) +
  geom_line(size = 1, color = "black") +
  geom_point(size = 2, color = "black") +
  scale_y_continuous(name = "AIC difference with Model 0") +
  scale_x_continuous(name = "Number of completed interviews", breaks = 1:20) +
  themeKUL + 
  theme(legend.position = "none",
        axis.line = element_line())

p

```


