---
title: "Estimation of Racial Disparities When Race is Not Observed"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Estimation of Racial Disparities When Race is Not Observed}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

This document will walk you through how to use `birdie`. First, load in the package.

```{r setup, message=FALSE}
library(birdie)
library(dplyr)
```

For a concrete example, we'll use some fake voter file data.
Our goal is to estimate turnout rates by race.

```{r}
data(pseudo_vf)

print(pseudo_vf)
```

You'll notice that we have a `race` column in our data.
That will allow us to check our work once we're done.
For now, we'll generate the *true* distribution of turnout given race, along with the marginal distribution of each variable.

```{r}
p_xr = prop.table(table(pseudo_vf$turnout, pseudo_vf$race), margin=2)
p_x = prop.table(table(pseudo_vf$turnout))
p_r = prop.table(table(pseudo_vf$race))
```

There are two steps to applying the `birdie` methodology:

1. Generate a first set of individual race probabilities using Bayesian Improved Surname Geocoding (`bisg()`).
1. For a specific outcome variable of interest, run a Bayesian Instrumental Regression for Disparity Estimation (`birdie()`) model to come up with estimated probabilities conditional on race.

## Generating BISG probabilities

For the first step, you can use any BISG software, including the [`wru` R package](https://cran.r-project.org/package=wru).
However, `birdie` provides its own `bisg()` function to make this easy and very computationally efficient.
To use `bisg()`, you provide a formula that labels the predictors used.
You use `nm()` to show which variable contains last names, which must always be provided.
ZIP codes and states can be labeled with `zip()` and `state()`.
Other types of geographies can be used as well---just read the documentation for `bisg()`.

```{r}
r_probs = bisg(~ nm(last_name) + zip(zip), data=pseudo_vf, p_r=p_r)
print(r_probs)
```

Each row `r_probs` matches a row in `pseudo_vf`.
**It's important to note that here we are assuming that we know the overall racial distribution of our population** (registered voters).
Because of that, we provide the `p_r=p_r` argument, which gives `bisg()` the overall racial distribution.
If you don't know the overall racial distribution in your context (even a guess is better than nothing), then you could pass in something like the national distribution of race (which is conveniently provided by `p_r_natl()`).

### Alternative race probabilities

Rather than predicting individual race with the standard BISG methodology, you may want to use the improved fully Bayesian Surname Improved Geocoding (fBISG) of [Imai et al. (2022)](https://doi.org/10.1126/sciadv.adc9824).
Compared to standard BISG, fBISG accounts for some of the measurement error in the Census counts.
This improves the calibration of the probabilities (which is important for accurate disparity estimation), and also improves accuracy among minority populations.

To generate fBISG probabilities, `birdie` provides the `bisg_me()` function, which works just like `bisg()`.

```{r}
r_probs_me = bisg_me(~ nm(last_name) + zip(zip), data=pseudo_vf, p_r=p_r, iter=2000)
```

Comparing to the standard BISG probabilities, the measurement-error-adjusted probabilities are often more calibrated.
One way to see this is to estimate the marginal distribution of race from the probabilities.

```{r}
colMeans(r_probs)
colMeans(r_probs_me)
# actual 
p_r
```

The fBISG probabilities are much closer to the actual distribution of race in the data than the standard BISG probabilities are.

### Why aren't BISG probabilities enough?

At this point, many analyses stop.
One can threshold the BISG probabilities to produce a single racial prediction for every individual.
Or one can use the BISG probabilities inside weighted averages and weighted regressions.

For example, we could try to estimate turnout rates by race, using the BISG probabilities as weights:

```{r}
est_weighted(r_probs, turnout ~ 1, data=pseudo_vf)
```

However, as discussed in the methodology paper ([McCartan et al. 2024](https://www.nber.org/papers/w32373)), **this approach is generally biased**.
Essentially, it only measures the part of the association between race and turnout that is mediated through names and locations.
It doesn't properly account for other ways in which race could be associated with the outcome.
The BIRDiE methodology addresses this problem by relying on **a different assumption: that names are independent of outcomes (here, turnout) conditional on location and race.**
For example, among White voters in a particular ZIP code, this assumption would mean that voters name Smith and those named Jones are both equally likely to vote.

## Estimating distributions by race

We're now ready to estimate turnout by race.
For this we'll use the `birdie()` function, and provide it with a formula describing our BIRDiE model, including our variable of interest `turnout` and our geography variable `zip`.
We provide `family=cat_mixed()` to indicate that we want to fit a Categorical mixed-effects regression model for turnout.
Here, we wrap `zip` in the `proc_zip()` function, which, among other things, recodes missing ZIP codes as "Other" so that the model doesn't encounter any missing data.
The first argument to `birdie()` is `r_probs`, the racial probabilities.
`birdie` knows how to handle its columns of because they came from this package.
If you use a different package, the columns may be named differently.
The `prefix` parameter to `birdie()` lets you specify the naming convention for your probabilities.

```{r messages=FALSE}
fit = birdie(r_probs, turnout ~ (1 | proc_zip(zip)), data=pseudo_vf, family=cat_mixed())
print(fit)
```

### Types of BIRDiE Models

The BIRDiE model we just fit is a *mixed-effects* model.
It estimates a different relationship between turnout and race in every ZIP, but partially pools these estimates towards a common global estimate of the turnout-race relationship.
**BIRDiE supports three other general types of models as well**: the complete-pooling and no-pooling categorical regression models, and a Normal linear model which can be used for continuous outcome variables (when the true regression function is assumed to be additive and linear in the covariates).
The complete-pooling model uses a formula like `turnout ~ 1` and only estimates a single, global relationship between turnout and race.
The model therefore assumes that turnout has no association with geography, after controlling for race.
The no-pooling model uses a formula like `turnout ~ proc_zip(zip)`.
While this model can be more computationally efficient to fit than the mixed-effects model, its performance can suffer on smaller datasets like the one used here.
We recommend the mixed-effects model for general use when the outcome variable is discrete.

### Extracting population and small-area  estimates

The `birdie()` function returns an object of class `birdie`, which supports many additional functions.
You can quickly extract the population turnout-race estimates using `coef()` or `tidy()`.
The former produces a matrix, while the latter returns a tidy data frame that may be useful in plotting or in downstream analyses.

```{r}
tidy(fit)
```

These estimates are quite close to the true distribution of turnout and race for most racial groups:

```{r}
coef(fit)
p_xr # Actual
```

The estimates suffer here for the smaller racial groups, which each comprise roughly 1-2% of the sample

You can also extract estimates by geography (and other covariates, if they are present in the model formula) by passing `subgroup=TRUE` to either `coef()` or `tidy()`.

```{r}
head(tidy(fit, subgroup=TRUE))
```

### Generating improved individual BISG probabilities

In addition to producing estimates for the whole sample and specific subgroups, BIRDiE yields improved individual race probabilities.
The "input" BISG probabilities are for race given surname and location.
The "output" probabilities from BIRDiE are for race given surname, location, and also turnout.
When the outcome variable is strongly associated with race, these BIRDiE-improved probabilities can be significantly more accurate than the standard BISG probabilities.

Accessing these improved probabilities is simple with the `fitted()` function.
```{r}
head(fitted(fit))

plot(r_probs$pr_white, fitted(fit)$pr_white, cex=0.1)
```

### Multiple Imputation from a BIRDiE Model

By simulating from the improved BISG probabilities, multiple imputations of the missing race assignments can be generated.
Each imputation can be fed through a downstream analysis, and the results combined by mixing posterior draws or via Rubin's rules.

For example, suppose we want to estimate the overall turnout rate for Black and Hispanic voters combined.
One way to do that would be as follows.
```{r}
race_lbl = levels(pseudo_vf$race)
calc_bh_turn = function(race_imp) {
    is_bh = race_lbl[race_imp] %in% c("black", "hisp")
    mean((pseudo_vf$turnout == "yes")[is_bh])
}

est_birdie = simulate(fit, 200) |> # 200 imputations stored as an integer matrix
    apply(2, calc_bh_turn) # calculate turnout for each imputation
hist(est_birdie)
```

We can compare the results with doing the same imputation procedure from the raw BISG probabilities.
The BIRDiE imputations are much closer to the true value.

```{r}
est_bisg = simulate(r_probs, 200) |> # simulate() works on BISG objects too
    apply(2, calc_bh_turn) 

tibble(
    actual = with(pseudo_vf, mean((turnout == "yes")[race %in% c("black", "hisp")])),
    est_birdie = mean(est_birdie),
    est_bisg = mean(est_bisg),
)
```

### Generating standard errors

One drawback of the computationally efficient EM algorithm that `birdie()` uses for model fitting is the lack of uncertainty quantification.
There are two approaches to generating standard errors for BIRDiE models, bootstrapping and Gibbs sampling, which are discussed below.
However, for most datasets, **non-sampling error in Census data and violations of model assumptions will cause much more bias than sampling variance**.

Gibbs sampling is preferred as it produces posterior draws from the full Bayesian model, but is not available for the mixed-effects model.
Calling `simulate()` to generate multiple imputations after Gibbs sampling will produce imputations that account for uncertainty in the model parameters, not just the missing data.
To use the Gibbs sampler for inference, provide `algorithm="gibbs"` to `birdie()`.
The `iter` parameter controls the number of bootstrap replicates.
Posterior variance estimates are accessible with `$se` or using the `vcov()` generic, and will be plotted with `plot()`.
The multiple imputation approach discussed above can also be used to understand uncertainty in other quantities of interest.

To bootstrap, simply set `algorithm="em_boot"` in `birdie()`.
Bootstrapping is also not yet available for the mixed-effects model.
The `iter` parameter controls the number of bootstrap replicates.
The standard errors are accessible with `$se` or using the `vcov()` generic,
and will be plotted with `plot()`.

```{r}
fit_boot = birdie(r_probs, turnout ~ 1, data=pseudo_vf, algorithm="em_boot", iter=200)

fit_boot$se
```