Intro to outcomerate

This vignette demonstrates the basic applications of the outcomerate package in R. I will draw on the popular tidyverse family of packages for the analysis.

# load packages
library(outcomerate)
library(dplyr)
library(tidyr)
library(knitr)

To keep things lighthearted, I will use a toy dataset named middleearth. The data consists of 1691 rows, each representing an attempt to interview a member of middle earth. Not all elements in the sample resulted in a completed interview, however. Some cases could not be located, others were located but no one was available, some individuals were found but refused to participate, etc. These particular ‘dispositions’ can be summarized from the code variable in the data:

# load dataset
data(middleearth)

# tabulate frequency table of outcomes
kable(count(middleearth, code, outcome))

It is common for survey practitioners to report a number of outcome rates. These rates give an indication as to the quality of the field work. For example, you may want to know the response rate: the proportion of all cases from our intended sample that actually resulted in an interview.

How might we go about calculating this?

When we inspect our disposition codes, it become apparent that there could be several ways to do this. For example, you may start by using the total number of complete cases (760) and diving this by the number of observations in the data, 760 / 1691 = 0.45. But what about partially completed interviews? If you include those, you would get a rate of (760 + 339) / 1691 = 0.65.

It turns out that there are a lot of ways to calculate such outcome rates. Unless we specify exactly what we mean by “response rate”, it is easy for claims regarding survey quality to become opaque, lacking comparability with other surveys. For this reason, the American Association for Public Opinion Research (AAPOR) has published a set of standardized definitions for practitioners. The guide has no fewer than 6 different variants of the ‘response rate.’ In the our example, the rates we calculated would match to AAPOR’s “Response Rate 1” and “Response Rate 2”:

$$ \textrm{RR1} = \frac{\textrm{I}}{\textrm{(I + P) + R + O + NC + (UO + UH)}} = \frac{760}{(760 + 339) + 59 + 1 + 288 + (173 + 0)} = 0.47 $$

$$ \textrm{RR2} = \frac{\textrm{(I + P)}}{\textrm{(I + P) + R + O + NC + (UO + UH)}} = \frac{(760 + 339)}{(760 + 339) + 59 + 1 + 288 + (173 + 0)} = 0.68 $$

What’s more, the guide has multiple definitions for contact rates, refusal rates, and cooperation rates, and weighted rates. It can easily become tedious to look all these up and calculate them by hand. The outcomerate package makes it easier by giving all rates (and more) in one go:

disp_counts <- c(I = 760, P = 339, R = 59, NC = 288, O = 1, UO = 173, NE = 71) 

e <- eligibility_rate(disp_counts)
outcomerate(disp_counts, e = e)
#>        RR1        RR2        RR3        RR4        RR5        RR6      COOP1 
#> 0.46913580 0.67839506 0.47149080 0.68180052 0.52522460 0.75950242 0.65573770 
#>      COOP2      COOP3      COOP4       REF1       REF2       REF3       CON1 
#> 0.94823123 0.65630397 0.94905009 0.03641975 0.03660258 0.04077402 0.71543210 
#>       CON2       CON3       LOC1       LOC2 
#> 0.71902347 0.80096752 0.89320988 0.89769367

Each of these rates has a precise definition (see ?outcomerate for details). As we can see, RR1 and RR2 match our earlier calculations. In the example, I needed to specify the parameter e, the estimated proportion of unknown cases unknowns (UO) that were eligible. The eligibility_rate() offers a default way to calculate this, but others may be appropriate.

If we had wanted just to return the two rates from above, we could specify this:

outcomerate(disp_counts, rate = c("RR1", "RR2"))
#>       RR1       RR2 
#> 0.4691358 0.6783951

# print the head of the dataset
head(middleearth)
#> # A tibble: 6 × 9
#>   code  outcome            researcher region    Q1       Q2   day race   svywt
#>   <ord> <ord>              <chr>      <fct>     <fct> <int> <int> <fct>  <dbl>
#> 1 UO    Unlocated          #23        Beleriand <NA>     NA     1 Elf       32
#> 2 I     Complete interview #23        Beleriand No        7     1 Hobbit    52
#> 3 I     Complete interview #23        Beleriand No        7     1 Hobbit    52
#> 4 P     Partial interview  #13        Beleriand No        7     1 Hobbit    52
#> 5 NE    Not eligible       #50        Beleriand <NA>     NA     1 Man       85
#> 6 I     Complete interview #23        Beleriand No        7     1 Man       85

# calculate rates using codes; should be same result as before
outcomerate(middleearth$code, e = e)
#>        RR1        RR2        RR3        RR4        RR5        RR6      COOP1 
#> 0.46913580 0.67839506 0.47149080 0.68180052 0.52522460 0.75950242 0.65573770 
#>      COOP2      COOP3      COOP4       REF1       REF2       REF3       CON1 
#> 0.94823123 0.65630397 0.94905009 0.03641975 0.03660258 0.04077402 0.71543210 
#>       CON2       CON3       LOC1       LOC2 
#> 0.71902347 0.80096752 0.89320988 0.89769367

# create a small wrapper function
get_rates <- function(x, ...){
  rlist <- c("RR1", "RR2", "COOP1", "COOP2", "CON1", "REF1", "LOC1")
  as.data.frame(as.list(outcomerate(x, rate = rlist, e = e, ...)))
}

# calculate rates by group
middleearth %>%
  group_by(race) %>%
  summarise(n     = n(),
            Nhat  = sum(svywt),
            rates = list(get_rates(code))) %>%
  unnest(cols = c(rates)) %>%
  kable(digits = 2, caption = "Outcome Rates by Race")

# calculate weighted rates by group
middleearth %>%
  group_by(region) %>%
  summarise(n     = n(),
            Nhat  = sum(svywt),
            rates = list(get_rates(code, weight = svywt))) %>%
  unnest(cols = c(rates)) %>%
  kable(digits = 2, caption = "Weighted Outcome Rates by Region")

library(ggplot2)
library(stringr)

# day-by-day quality monitoring
middleearth %>%
  group_by(day) %>%
  summarise(rates = list(get_rates(code))) %>%
  unnest(cols = c(rates)) %>%
  gather(rate, value, -day) %>%
  mutate(type = str_sub(rate, start = -9, end = -2)) %>%
  ggplot(aes(x = day, y = value, colour = rate)) +
  geom_line(size = 1) +
  facet_wrap(~type) +
  labs(title = "Outcome Rates Over Time")
#> Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
#> ℹ Please use `linewidth` instead.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
#> generated.

# first, calculate the outcome rates
(res <- outcomerate(middleearth$code))
#>        RR1        RR2        RR5        RR6      COOP1      COOP2      COOP3 
#> 0.46913580 0.67839506 0.52522460 0.75950242 0.65573770 0.94823123 0.65630397 
#>      COOP4       REF1       REF3       CON1       CON3       LOC1 
#> 0.94905009 0.03641975 0.04077402 0.71543210 0.80096752 0.89320988

# estimate standard errors using the Normal approximation for proportions 
se <- sapply(res, function(p) sqrt((p * (1 - p)) / nrow(middleearth)))

# calculate 95% confidence intervals
rbind(res - (se * 1.96), res + (se * 1.96))
#>            RR1       RR2       RR5       RR6     COOP1     COOP2     COOP3
#> [1,] 0.4453496 0.6561319 0.5014233 0.7391318 0.6330916 0.9376710 0.6336667
#> [2,] 0.4929220 0.7006582 0.5490259 0.7798730 0.6783838 0.9587915 0.6789412
#>          COOP4       REF1       REF3      CON1      CON3      LOC1
#> [1,] 0.9385691 0.02749088 0.03134782 0.6939260 0.7819369 0.8784892
#> [2,] 0.9595310 0.04534863 0.05020021 0.7369382 0.8199982 0.9079305