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

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.

data(pseudo_vf)

print(pseudo_vf)
#> # A tibble: 5,000 × 4
#>    last_name zip   race  turnout
#>    <fct>     <fct> <fct> <fct>
#>  1 BEAVER    28748 white yes
#>  2 WILLIAMS  28144 black no
#>  3 ROSEN     28270 white yes
#>  4 SMITH     28677 black yes
#>  5 FAY       28748 white no
#>  6 CHURCH    28215 white yes
#>  7 JOHNSON   28212 black yes
#>  8 SZCZYGIEL NA    white yes
#>  9 SUMMERS   28152 black yes
#> 10 STARLING  28650 white yes
#> # … with 4,990 more rows

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.

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()).
2. 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. 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_probs = bisg(~ nm(last_name) + zip(zip), data=pseudo_vf, p_r=p_r)
print(r_probs)
#> # A tibble: 5,000 × 6
#>    pr_white pr_black  pr_hisp pr_asian  pr_aian pr_other
#>       <dbl>    <dbl>    <dbl>    <dbl>    <dbl>    <dbl>
#>  1    0.971  0.00517 0.00201  0.000144 0.0136    0.00819
#>  2    0.128  0.860   0.00186  0.000169 0.00104   0.00920
#>  3    0.979  0.00539 0.00436  0.00232  0.000606  0.00837
#>  4    0.521  0.459   0.00535  0.000219 0.00150   0.0124
#>  5    0.989  0.00165 0.00257  0.000319 0.00182   0.00506
#>  6    0.522  0.431   0.0175   0.00125  0.00384   0.0242
#>  7    0.112  0.859   0.0113   0.000720 0.00190   0.0153
#>  8    0.991  0.00794 0.000766 0        0         0
#>  9    0.701  0.283   0.00183  0.000160 0.00136   0.0122
#> 10    0.854  0.133   0.00265  0.000204 0.00241   0.00845
#> # … with 4,990 more rows

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). Compared to standard BISG, fBISG accounts for 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_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 measurment-error-adjusted probabilities are often more calibrated. One way to see this is to estimate the marginal distribution of race from the probabilities.

colMeans(r_probs)
#>    pr_white    pr_black     pr_hisp    pr_asian     pr_aian    pr_other
#> 0.656372686 0.277807934 0.032031532 0.010139008 0.009564482 0.014084358
colMeans(r_probs_me)
#>  pr_white  pr_black   pr_hisp  pr_asian   pr_aian  pr_other
#> 0.7089420 0.2001877 0.0433253 0.0186468 0.0118867 0.0170115
# actual
p_r
#>
#>  white  black   hisp  asian   aian  other
#> 0.7178 0.2078 0.0338 0.0114 0.0102 0.0190

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:

est_weighted(r_probs, turnout ~ 1, data=pseudo_vf)
#> Weighted estimator
#> Formula: turnout ~ 1
#>    Data: pseudo_vf
#> Number of obs: 5,000; groups: 1
#> Estimated distribution:
#>     white black  hisp asian  aian other
#> no  0.316 0.335 0.385 0.531 0.502 0.341
#> yes 0.684 0.665 0.615 0.469 0.498 0.659

However, as discussed in the methodology paper (McCartan et al. 2023), this approach is generally biased. Essentially, it only measures the association between race and turnout which 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. 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.

fit = birdie(r_probs, turnout ~ (1 | proc_zip(zip)), data=pseudo_vf)
#> Using default prior for Pr(X | R):
#> → Prior scale on fixed effects coefficients: 1.0
#> → Prior mean of random effects standard deviation: 0.20
#> This message is displayed once per session.
print(fit)
#> Multinomial mixed-effects BIRDiE model
#> Formula: turnout ~ (1 | proc_zip(zip))
#>    Data: pseudo_vf
#> Number of obs: 5,000; groups: 626
#> Estimated distribution:
#>     white black  hisp asian  aian other
#> no  0.299 0.363 0.418 0.652 0.698 0.228
#> yes 0.701 0.637 0.582 0.348 0.302 0.772

### 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 two other general types of models as well: the complete-pooling and no-pooling models. The former 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 latter 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 perfromance can suffer on smaller datasets like the one used here. We recommend the mixed-effects model for general use.

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

tidy(fit)
#> # A tibble: 12 × 3
#>    turnout race  estimate
#>    <chr>   <chr>    <dbl>
#>  1 no      white    0.299
#>  2 yes     white    0.701
#>  3 no      black    0.363
#>  4 yes     black    0.637
#>  5 no      hisp     0.418
#>  6 yes     hisp     0.582
#>  7 no      asian    0.652
#>  8 yes     asian    0.348
#>  9 no      aian     0.698
#> 10 yes     aian     0.302
#> 11 no      other    0.228
#> 12 yes     other    0.772

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

coef(fit)
#>         white     black      hisp     asian      aian     other
#> no  0.2993663 0.3630225 0.4180173 0.6515683 0.6975264 0.2282492
#> yes 0.7006337 0.6369775 0.5819827 0.3484317 0.3024736 0.7717508
p_xr # Actual
#>
#>           white     black      hisp     asian      aian     other
#>   no  0.3014767 0.3570741 0.4142012 0.4385965 0.6666667 0.5894737
#>   yes 0.6985233 0.6429259 0.5857988 0.5614035 0.3333333 0.4105263

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

head(tidy(fit, subgroup=TRUE))
#> # A tibble: 6 × 4
#>   zip   turnout race  estimate
#>   <fct> <chr>   <chr>    <dbl>
#> 1 28748 no      white    0.269
#> 2 28748 yes     white    0.731
#> 3 28748 no      black    0.373
#> 4 28748 yes     black    0.627
#> 5 28748 no      hisp     0.411
#> 6 28748 yes     hisp     0.589

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

head(fitted(fit))
#> # A tibble: 6 × 6
#>   pr_white pr_black pr_hisp  pr_asian  pr_aian pr_other
#>      <dbl>    <dbl>   <dbl>     <dbl>    <dbl>    <dbl>
#> 1  0.979    0.00447 0.00164 0.0000702 0.00584   0.00862
#> 2  0.00816  0.986   0.00131 0.000182  0.00120   0.00364
#> 3  0.971    0.00920 0.00480 0.00166   0.000382  0.0129
#> 4  0.525    0.457   0.00413 0.000113  0.000680  0.0138
#> 5  0.984    0.00228 0.00390 0.000761  0.00464   0.00444
#> 6  0.562    0.399   0.0160  0.000495  0.00134   0.0207

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

### Bootstrapping standard errors

One drawback of the computationally efficient EM algorithm that birdie() uses for model fitting is the lack of uncertainty quantification. The function does however support bootstrapping for some model types. In general, we recommend bootstrapping a simple complete-pooling model to understand the rough amount of sampling variation present. For most datasets, non-sampling error in Census data and violations of model assumptions will cause much more bias than sampling variance.

To bootstrap, simply pass a nonzero value of se_boot to birdie(). This parameter controls the number of bootstrap replicates. The standard errors are accessible with $se or using the vcov() generic. fit_boot = birdie(r_probs, turnout ~ 1, data=pseudo_vf, se_boot=200) #> Using c(1+ε, 1+ε, ..., 1+ε) prior for Pr(X | R) #> This message is displayed once per session. fit_boot$se
#>         white      black       hisp      asian      aian     other
#> no  0.0105847 0.02161647 0.05131842 0.09558014 0.1124541 0.3518539
#> yes 0.0105847 0.02161647 0.05131842 0.09558014 0.1124541 0.3518539