Observing interviewer performance in slices or by traces: A comparison of methods to predict interviewers’ individual contributions to interviewer variance

Celine Wuyts and Geert Loosveldt

26 februari 2022

1 Data and methods

# 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)

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))

1.1 Selection of survey 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))))

##        
##         FALSE TRUE Sum
##   FALSE     0   42  42
##   TRUE     71   68 139
##   Sum      71  110 181

# 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))))

##        
##         FALSE TRUE Sum
##   FALSE     0   22  22
##   TRUE     36    3  39
##   Sum      36   25  61

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))

Intra-interviewer correlations for selected survey items
		Intra-interviewer correlation
Name	Label	ESS6 BEDUT	ESS7 BEDUT
iprspot	Important to get respect from others	0.104	0.097
prtdgcl	How close to party	0.087	.
dfegcon	Different race or ethnic group: contact, how often	.	0.086
tvpol	TV watching, news, politics, current affairs on average weekday	.	0.085
dspplvt	Voters discuss politics with people they know before deciding how to vote	0.082	.
iplylfr	Important to be loyal to friends and devote to people close	0.076	.
freehms	Gays and lesbians free to live life as they wish	0.073	.
gvctzpv	The government protects all citizens against poverty	0.073	.
qfimedu	Qualification for immigration: good educational qualifications	.	0.069
cttresa	The courts treat everyone the same	0.068	.
gvexpdc	The government explains its decisions to voters	0.067	.
votedir	Citizens have the final say on political issues by voting directly in referendums	0.067	.
alcbnge	Frequency of binge drinking for men and women, last 12 months	.	0.067
gvcodmc	In country government formed by coalition	0.064	.
ipsuces	Important to be successful and that people recognize achievements	0.060	.
euftf	European Union: European unification go further or gone too far	.	0.059
ipshabt	Important to show abilities and be admired	0.059	.
dfprtal	Different political parties offer clear alternatives to one another	0.054	.
impenv	Important to care for nature and environment	0.053	.
qfimlng	Qualification for immigration: speak country’s official language	.	0.052
gptpelcc	In country governing parties are punished in elections when they have done a bad job	0.052	.
flclpla	Feel close to the people in local area	0.052	.
tmendng	Enthusiastic about what you are doing, how much of the time	0.052	.
tnapsur	Take notice of and appreciate your surroundings	0.052	.
eimpcnt	Allow many, few immigrants from poorer countries in Europe	.	0.051
grdfincc	In country the government takes measures to reduce differences in income levels	0.050	.
imbleco	Taxes and services: immigrants take out more than they put in or less	.	0.049
slprl	Sleep was restless, how often past week	0.048	.
algyplv	Allow many or few Gypsies to come and live in country	.	0.048
noimbro	Of every 100 people in country how many born outside country	.	0.047
imdfetn	Allow many, few immigrants of different race, ethnic group from majority	.	0.044
ipstrgv	Important that government is strong and ensures safety	0.039	0.049
gvctzpvc	In country the government protects all citizens against poverty	0.043	.
qfimchr	Qualification for immigration: Christian background	.	0.043
impdiff	Important to try new and different things in life	.	0.042
plinsoc	Your place in society	0.042	.
rghmgpr	The rights of minority groups are protected	0.041	.
ctstogv	The courts able to stop the government acting beyond its authority	0.041	.
impfun	Important to seek fun and things that give pleasure	0.040	0.041
gvrfgap	Government should be generous judging applications for refugee status	.	0.040
deaimpp	Deal with important problems in life	0.040	.
etapapl	Easy to take part in politics	.	0.040
cptppol	Confident in own ability to participate in politics	.	0.039
imwbcnt	Immigrants make country worse or better place to live	.	0.039
pltaviec	In country politicians take into account the views of other European governments	0.039	.
trtrsp	Feel people treat you with respect	0.038	.
gincdif	Government should reduce differences in income levels	0.038	.
lotsgot	There are lots of things I am good at	0.038	.
votedirc	In country citizens have the final say on political issues by voting directly in referendums	0.037	.
ipmodst	Important to be humble and modest, not draw attention	0.037	.
ipeqopt	Important that people are treated equally and have equal opportunities	.	0.036
prhlppl	Provide help and support to people you are close to	0.036	.
ipgdtim	Important to have a good time	0.036	.
ipadvnt	Important to seek adventures and have an exciting life	.	0.035
ipudrst	Important to understand different people	0.035	.
pltavie	Politicians take into account the views of other European governments	0.035	.
happy	How happy are you	0.034	.
psppsgv	Political system allows people to have a say in what government does	.	0.034
enjlf	Enjoyed life, how often past week	.	0.032
flapppl	Feel appreciated by people you are close to	0.032	.
gvtrimg	Compared to yourself government treats new immigrants better or worse	.	0.031

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

**Descriptives of intra-interviewer correlations (all survey items)**

Statistic	N	Mean	St. Dev.	Min	Pctl(25)	Median	Pctl(75)	Max

icc_fit_controls	249	0.023	0.021	0.000	0.006	0.020	0.035	0.104

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)")

**Descriptives of intra-interviewer correlations (all survey items)**

Statistic	N	Mean	St. Dev.	Min	Pctl(25)	Median	Pctl(75)	Max

ESS6_BEDUT	139	0.025	0.022	0.000	0.007	0.023	0.037	0.104
ESS7_BEDUT	110	0.021	0.020	0.000	0.005	0.017	0.031	0.097

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

**Descriptives of intra-interviewer correlations (selected survey items)**

Statistic	N	Mean	St. Dev.	Min	Pctl(25)	Median	Pctl(75)	Max

icc_fit_controls	64	0.051	0.017	0.031	0.039	0.046	0.059	0.104

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)")

**Descriptives of intra-interviewer correlations (selected survey items)**

Statistic	N	Mean	St. Dev.	Min	Pctl(25)	Median	Pctl(75)	Max

ESS6_BEDUT	39	0.052	0.017	0.032	0.038	0.048	0.062	0.104
ESS7_BEDUT	25	0.050	0.018	0.031	0.039	0.044	0.052	0.097

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)

1.2 Interviewers’ deviation from the standardized interviewing protocol

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)

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)

Checklist of standardized interviewing protocol
	Deviation from check
Check	ESS6 BEDUT	ESS7 BEDUT
Reads introduction	12.5%	11.2%
Speaks the language of the interview fluently	0.0%	7.5%
Reads questions completely	36.4%	51.2%
Does not add anything to questions	15.9%	32.5%
Reads all applicable questions	4.5%	8.8%
Reads introductory sentences before questions	38.6%	27.5%
Repeats questions in case of irrelevant or unclear answers or when requested by R	9.1%	20.0%
Does not read interviewer instructions	6.8%	26.2%
Reads references to showcards	5.7%	15.0%
Does not read showcards	38.6%	33.8%
Asks for additional explanation in case answer is not one of the available options	36.4%	28.7%
Does not provide example answers	27.3%	17.5%
Does not read Refusal, Dont know, and Other	1.1%	1.2%
Does not read additional options within brackets	2.3%	2.5%
Reads all options in case no show card available	20.5%	33.8%
Probes at least once in case of Refusal or Dont know	4.5%	2.5%
Probes at least once in case off all-that-apply questions	30.7%	15.0%
Is not suggestive or steering	29.5%	22.5%
Does not give his/her opinion	3.4%	3.8%
Asks R to interpret question him/herself in case R asks for explanation	6.8%	21.2%
Does not read questions too slow or too fast	1.1%	16.2%
Does not read questions too loud or too quiet	0.0%	1.2%
Reads questions clearly	4.5%	6.2%
Is agreeable to listen to	3.4%	6.2%
Reads all questions in the same way, without apology	5.7%	5.0%
Gives R sufficient time to answer	5.7%	6.2%
Gives short confirmations	0.0%	0.0%
Is friendly and interested	0.0%	2.5%
Does not give value judgements, approvals, disapprovals	1.1%	1.2%

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

**Descriptives of standardized interviewing deviation score**

Statistic	N	Mean	St. Dev.	Min	Pctl(25)	Pctl(75)	Max

ESS6_BEDUT	88	3.523	3.051	0.000	1.000	5.000	15.000
ESS7_BEDUT	80	4.275	3.670	0.000	2.000	5.250	16.000

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

##           essround
## AUDIO_NDEV          6          7
##         0  0.12500000 0.08750000
##         1  0.18181818 0.12500000
##         2  0.14772727 0.13750000
##         3  0.14772727 0.21250000
##         4  0.07954545 0.10000000
##         5  0.10227273 0.08750000
##         6  0.03409091 0.03750000
##         7  0.09090909 0.01250000
##         8  0.02272727 0.07500000
##         9  0.02272727 0.03750000
##         10 0.01136364 0.00000000
##         11 0.01136364 0.02500000
##         12 0.01136364 0.01250000
##         13 0.00000000 0.01250000
##         14 0.00000000 0.01250000
##         15 0.01136364 0.01250000
##         16 0.00000000 0.01250000

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

##        essround
##             6     7
##   FALSE 0.500 0.275
##   TRUE  0.500 0.725

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

##        essround
##                 6         7
##   FALSE 0.4545455 0.4625000
##   TRUE  0.5454545 0.5375000

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

##        essround
##                  6          7
##   FALSE 0.93181818 0.82500000
##   TRUE  0.06818182 0.17500000

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

Distribution of the interviewers’ standardized interviewing deviation score

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

Joint distribution of the two survey rounds’ standardized interviewing deviation score

1.3 Interviewers’ interview speed

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

**Descriptives of average interview speed**

Statistic	N	Mean	St. Dev.	Min	Pctl(25)	Pctl(75)	Max

ESS6_BEDUT	88	4.197	0.481	2.398	3.914	4.522	5.058
ESS7_BEDUT	80	3.588	0.527	2.075	3.290	3.855	4.866

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

**Descriptives of average interview acceleration over the questionnaire**

Statistic	N	Mean	St. Dev.	Min	Pctl(25)	Pctl(75)	Max

ESS6_BEDUT	88	0.103	0.158	-0.296	0.003	0.171	0.654
ESS7_BEDUT	80	-0.150	0.140	-0.503	-0.249	-0.055	0.224

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

Joint distribution of the two survey rounds’ average interview 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

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

1.4 Interviewers’ contributions to measurement error

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"))

# Deviations larger than 0.5 in absolute value

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

## 
## FALSE  TRUE 
##  5411    19

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

## 
##       FALSE        TRUE 
## 0.996500921 0.003499079

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

Interviewers’ individual absolute bias estimates

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

2 Results

# 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))

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)))

**Fixed parameter estimates for models explaining interviewers’ contributions to interviewer variance**

	Model 0	Model A	Model A’	Model B	Model C	Model C’
Standardized interviewing deviation count		0.0547^***			0.0553^***
		(0.0202)			(0.0201)

At least one flag for question reading			-0.0310			-0.0371^*
			(0.0204)			(0.0205)

At least one flag for deviation from neutrality			0.0463^**			0.0511^**
			(0.0203)			(0.0201)

At least one flag for inadequate interview pace			0.0460^**			0.0456^**
			(0.0193)			(0.0191)

Average interview speed				0.0031	0.0138	-0.0001
				(0.0250)	(0.0253)	(0.0259)

Average acceleration over questionnaire				0.1073^***	0.1075^***	0.1130^***
				(0.0243)	(0.0244)	(0.0247)

Number of interviews completed by interviewer	0.1691^***	0.1773^***	0.1612^***	0.1612^***	0.1678^***	0.1527^***
	(0.0201)	(0.0203)	(0.0209)	(0.0202)	(0.0204)	(0.0209)

Round 7	-0.0059	-0.0088	-0.0015	0.0622^*	0.0648^*	0.0702^**
	(0.0280)	(0.0280)	(0.0287)	(0.0346)	(0.0346)	(0.0354)

Constant	-0.0126	-0.0161	-0.0145	-0.0135	-0.0171	-0.0156
	(0.0418)	(0.0418)	(0.0422)	(0.0413)	(0.0414)	(0.0417)


Interviewer variance	0.055600	0.055006	0.059217	0.050811	0.051425	0.054879
Item variance	0.068541	0.068546	0.068560	0.068552	0.068561	0.068578
Residual variance	0.833192	0.832340	0.830945	0.831723	0.830527	0.829054
AIC	14,720	14,714	14,714	14,704	14,699	14,697
Observations	5,430	5,430	5,430	5,430	5,430	5,430

Note:	p<0.1; p<0.05; p<0.01

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))

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