Intro to outcomerate
In outcomerate: AAPOR Survey Outcome Rates

knitr::opts_chunk$set(
  fig.width = 7,
  fig.height = 4,
  collapse = TRUE,
  comment = "#>"
)

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

attach(as.list(table(middleearth$code)))

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 (r I) and diving this by the number of observations in the data, r I / r nrow(middleearth) = r round(I / nrow(middleearth), 2). But what about partially completed interviews? If you include those, you would get a rate of (r I + r P) / r nrow(middleearth) = r round((I + P) / nrow(middleearth), 2).

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{r I}{(r I + r P) + r R + r O + r NC + (r UO + 0)} = r round(outcomerate(middleearth$code, rate = "RR1"), 2) $$

$$ \textrm{RR2} = \frac{\textrm{(I + P)}}{\textrm{(I + P) + R + O + NC + (UO + UH)}} = \frac{(r I + r P)}{(r I + r P) + r R + r O + r NC + (r UO + 0)} = r round(outcomerate(middleearth$code, rate = "RR2"), 2) $$

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)

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

More Advanced Uses

In certain situations, you may want to calculate outcome rates based on a vector of codes, rather than a table of frequency counts. It is just as easy to obtain rates this way using outcomerate:

# print the head of the dataset
head(middleearth)

# calculate rates using codes; should be same result as before
outcomerate(middleearth$code, e = e)

Why might we prefer this input format, when it is just as easy to specify the counts?

Well, if we want to calculate outcome rates by some other covariate, we typically need to go back to the original data. For example, here we use dplyr and tidyr to calculate outcome rates of interest by race:

# 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() %>%
  kable(digits = 2, caption = "Outcome Rates by Race")

Weighted Outcome Rates

In certain situations, we also wish to produce weighted outcome rates, using the survey weights that are provided in the data. This is easy to do with one additional parameter:

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

Compare this to the equivalent unweighted estimates, and you see that the results are not the same.

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

By Date

Lastly, another useful application of grouped analysis is to calculate the rates by date. This allows you to monitor the quality day by day and notice if performance starts to change over time.

library(ggplot2)
library(stringr)

# day-by-day quality monitoring
middleearth %>%
  group_by(day) %>%
  summarise(rates = list(get_rates(code))) %>%
  unnest() %>%
  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")

In this example, we can see that the contact rate (CON) and response rate (RR) start to degrade in quality towards day 30. If fieldwork was still continuing, this could be something to look into and attempt to explain and/or redress.

Variance Estimation

To estimate the errors from estimates generated by outcomerate(), the simplest approach is to use the normal appromation. Since outcome rates are nothing more than proportions (or nearly so), their standard error is given by $SE(p) = \sqrt{(p(1-p))/n}$.

# first, calculate the outcome rates
(res <- outcomerate(middleearth$code))

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

With the standard error in hand, we can then construct frequentist confidence intervals:

# calculate 95% confidence intervals
rbind(res - (se * 1.96), res + (se * 1.96))

Weighted variance estimation in complex surveys require different procedures that go beyond the scope of this vignette. We recommend using svycontrast() from the survey package to obtain design-based errors that account for elements such as clustering and stratification. Bootstrapping primary sampling units (PSUs) may also be an appropriate method depending on the design at hand.