Example 09: Databases and APIs

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.2.1     ✔ readr     2.2.0
✔ forcats   1.0.1     ✔ stringr   1.6.0
✔ ggplot2   4.0.3     ✔ tibble    3.3.1
✔ lubridate 1.9.5     ✔ tidyr     1.3.2
✔ purrr     1.2.2     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
here() starts at /Users/kjhealy/Documents/courses/mptc

Working with
tidycensus

Load our packages

library(here)       # manage file paths
library(socviz)     # data and some useful things, especially %notin%
library(tidyverse)  # your friend and mine

library(scales)     # Convenient scale labels

Attaching package: 'scales'
The following object is masked from 'package:purrr':

    discard
The following object is masked from 'package:readr':

    col_factor
library(tidycensus) # Tidily talk to the Census
library(sf)         # Draw maps with ggplot
Linking to GEOS 3.13.0, GDAL 3.8.5, PROJ 9.5.1; sf_use_s2() is TRUE
# Don't needlessly download geometries files multiple times
options(tigris_use_cache = TRUE)

Census Data

  • Make sure you have a Census API key. Load it after you load tidycensus
library(tidycensus)

census_api_key("1234567890987654321")

Decennial Census Variables

There are a lot of them:

census_vars <- load_variables(year = 2010,
                              dataset = "sf1",
                              cache = TRUE)
Warning: `cache` is deprecated and ignored. tidycensus no longer writes
variable metadata to a local cache.
census_vars
# A tibble: 8,959 × 3
   name    label                                concept         
   <chr>   <chr>                                <chr>           
 1 H001001 Total                                HOUSING UNITS   
 2 H002001 Total                                URBAN AND RURAL 
 3 H002002 Total!!Urban                         URBAN AND RURAL 
 4 H002003 Total!!Urban!!Inside urbanized areas URBAN AND RURAL 
 5 H002004 Total!!Urban!!Inside urban clusters  URBAN AND RURAL 
 6 H002005 Total!!Rural                         URBAN AND RURAL 
 7 H002006 Total!!Not defined for this file     URBAN AND RURAL 
 8 H003001 Total                                OCCUPANCY STATUS
 9 H003002 Total!!Occupied                      OCCUPANCY STATUS
10 H003003 Total!!Vacant                        OCCUPANCY STATUS
# ℹ 8,949 more rows

Median Age by State:

age10 <- get_decennial(geography = "state",
                       variables = "P013001",
                       year = 2010)
Getting data from the 2010 decennial Census
Using Census Summary File 1
age10
# A tibble: 52 × 4
   GEOID NAME        variable value
   <chr> <chr>       <chr>    <dbl>
 1 01    Alabama     P013001   37.9
 2 02    Alaska      P013001   33.8
 3 04    Arizona     P013001   35.9
 4 05    Arkansas    P013001   37.4
 5 06    California  P013001   35.2
 6 22    Louisiana   P013001   35.8
 7 21    Kentucky    P013001   38.1
 8 08    Colorado    P013001   36.1
 9 09    Connecticut P013001   40  
10 10    Delaware    P013001   38.8
# ℹ 42 more rows

Some County Data for North Carolina

# Census variable names
popvars <- c("P005003", "P005004", "P005006", "P004003")


# Get a county-level dataset for NC with these variables
# The summary value is the total population.
nc <- get_decennial(geography = "county",
                    variables = popvars,
                    year = 2010,
                    summary_var = "P001001",
                    state = "NC") |>
  mutate(pct = 100 * (value / summary_value))
Getting data from the 2010 decennial Census
Using Census Summary File 1
nc
# A tibble: 400 × 6
   GEOID NAME                             variable  value summary_value   pct
   <chr> <chr>                            <chr>     <dbl>         <dbl> <dbl>
 1 37007 Anson County, North Carolina     P005003   12344         26948  45.8
 2 37011 Avery County, North Carolina     P005003   16029         17797  90.1
 3 37003 Alexander County, North Carolina P005003   32671         37198  87.8
 4 37015 Bertie County, North Carolina    P005003    7393         21282  34.7
 5 37013 Beaufort County, North Carolina  P005003   31705         47759  66.4
 6 37005 Alleghany County, North Carolina P005003    9862         11155  88.4
 7 37001 Alamance County, North Carolina  P005003  101718        151131  67.3
 8 37009 Ashe County, North Carolina      P005003   25420         27281  93.2
 9 37017 Bladen County, North Carolina    P005003   19242         35190  54.7
10 37019 Brunswick County, North Carolina P005003   86818        107431  80.8
# ℹ 390 more rows

We can do this a little more elegantly:

# Census variable names
popvars <- tribble(
  ~variable, ~name,
  "P005003", "nh_white",
  "P005004", "nh_black",
  "P005006", "nh_asian",
  "P004003", "hispanic")

popvars
# A tibble: 4 × 2
  variable name    
  <chr>    <chr>   
1 P005003  nh_white
2 P005004  nh_black
3 P005006  nh_asian
4 P004003  hispanic
# Get a county-level dataset for NC with these variables
# The summary value is the total population.
nc <- get_decennial(geography = "county",
                    variables = popvars$variable, 
                    year = 2010,
                    summary_var = "P001001",
                    state = "NC") |>
  mutate(pct = 100 * (value / summary_value))
Getting data from the 2010 decennial Census
Using Census Summary File 1
nc
# A tibble: 400 × 6
   GEOID NAME                             variable  value summary_value   pct
   <chr> <chr>                            <chr>     <dbl>         <dbl> <dbl>
 1 37007 Anson County, North Carolina     P005003   12344         26948  45.8
 2 37011 Avery County, North Carolina     P005003   16029         17797  90.1
 3 37003 Alexander County, North Carolina P005003   32671         37198  87.8
 4 37015 Bertie County, North Carolina    P005003    7393         21282  34.7
 5 37013 Beaufort County, North Carolina  P005003   31705         47759  66.4
 6 37005 Alleghany County, North Carolina P005003    9862         11155  88.4
 7 37001 Alamance County, North Carolina  P005003  101718        151131  67.3
 8 37009 Ashe County, North Carolina      P005003   25420         27281  93.2
 9 37017 Bladen County, North Carolina    P005003   19242         35190  54.7
10 37019 Brunswick County, North Carolina P005003   86818        107431  80.8
# ℹ 390 more rows
nc |>
  left_join(popvars, by = "variable") |>
  relocate(name, .after = variable) |>
  select(-variable) |>
  pivot_wider(names_from = name,
              names_glue = "{name}_{.value}",
              values_from = c(value, pct)) |>
  mutate(NAME = str_remove(NAME, " County.*"))
# A tibble: 100 × 11
   GEOID NAME      summary_value nh_white_value nh_black_value nh_asian_value
   <chr> <chr>             <dbl>          <dbl>          <dbl>          <dbl>
 1 37007 Anson             26948          12344          13038            281
 2 37011 Avery             17797          16029            706             56
 3 37003 Alexander         37198          32671           2025            357
 4 37015 Bertie            21282           7393          13252            103
 5 37013 Beaufort          47759          31705          12105            158
 6 37005 Alleghany         11155           9862            133             53
 7 37001 Alamance         151131         101718          27985           1806
 8 37009 Ashe              27281          25420            148            105
 9 37017 Bladen            35190          19242          12202             67
10 37019 Brunswick        107431          86818          12120            560
# ℹ 90 more rows
# ℹ 5 more variables: hispanic_value <dbl>, nh_white_pct <dbl>,
#   nh_black_pct <dbl>, nh_asian_pct <dbl>, hispanic_pct <dbl>

Put it in an object:

Now we have:

nc_wide
# A tibble: 100 × 11
   GEOID NAME      summary_value nh_white_value nh_black_value nh_asian_value
   <chr> <chr>             <dbl>          <dbl>          <dbl>          <dbl>
 1 37007 Anson             26948          12344          13038            281
 2 37011 Avery             17797          16029            706             56
 3 37003 Alexander         37198          32671           2025            357
 4 37015 Bertie            21282           7393          13252            103
 5 37013 Beaufort          47759          31705          12105            158
 6 37005 Alleghany         11155           9862            133             53
 7 37001 Alamance         151131         101718          27985           1806
 8 37009 Ashe              27281          25420            148            105
 9 37017 Bladen            35190          19242          12202             67
10 37019 Brunswick        107431          86818          12120            560
# ℹ 90 more rows
# ℹ 5 more variables: hispanic_value <dbl>, nh_white_pct <dbl>,
#   nh_black_pct <dbl>, nh_asian_pct <dbl>, hispanic_pct <dbl>

Top 20 Counties by Percent Hispanic

nc_wide |>
  slice_max(order_by = hispanic_pct, n = 20) |>
  ggplot(mapping = aes(x = hispanic_pct,
                       y = reorder(NAME, hispanic_pct))) +
  geom_point() +
  labs(x = "Percent Hispanic",
       y = NULL,
       title = "North Carolina Counties")

Log population vs Percent Hispanic

nc_wide |>
  ggplot(mapping = aes(x = summary_value,
                       y = hispanic_pct)) +
  geom_point() +
  scale_x_log10(labels = label_number(scale_cut = cut_short_scale())) +
  labs(x = "Log Population",
       y = "Percent Hispanic",
       title = "North Carolina Counties")

Log population vs Percent Black

nc_wide |>
  ggplot(mapping = aes(x = summary_value,
                       y = nh_black_pct)) +
  geom_point() +
  scale_x_log10(labels = label_number(scale_cut = cut_short_scale())) +
  labs(x = "Log Population",
       y = "Percent Black",
       title = "North Carolina Counties")

Percent White vs Percent Black

nc_wide |>
  ggplot(mapping = aes(x = nh_white_pct,
                       y = nh_black_pct)) +
  geom_point() +
  labs(x = "Percent White",
       y = "Percent Black",
       title = "North Carolina Counties")

Population Pyramids

We can ask the Census for its estimates of the age breakdown of the population, using the get_estimates() function.

usa <- get_estimates(
    geography = "us",
    product = "characteristics",
    breakdown = c("SEX", "AGEGROUP"),
    breakdown_labels = TRUE,
    year = 2019
  )

usa
# A tibble: 96 × 5
   GEOID NAME              value SEX        AGEGROUP          
   <chr> <chr>             <dbl> <chr>      <fct>             
 1 1     United States 328239523 Both sexes All ages          
 2 1     United States  19576683 Both sexes Age 0 to 4 years  
 3 1     United States  20195895 Both sexes Age 5 to 9 years  
 4 1     United States 161657324 Male       All ages          
 5 1     United States  10009207 Male       Age 0 to 4 years  
 6 1     United States  10322762 Male       Age 5 to 9 years  
 7 1     United States 166582199 Female     All ages          
 8 1     United States   9567476 Female     Age 0 to 4 years  
 9 1     United States   7528626 Female     Age 70 to 74 years
10 1     United States  20798268 Both sexes Age 10 to 14 years
# ℹ 86 more rows

Population Pyramids

usa_pyr <- usa |>
  filter(AGEGROUP %notin% "All ages",
         str_detect(AGEGROUP, "Age"),
         SEX %notin% "Both sexes") |>
  mutate(value = ifelse(SEX == "Male", -value, value))

usa_pyr
# A tibble: 36 × 5
   GEOID NAME              value SEX    AGEGROUP          
   <chr> <chr>             <dbl> <chr>  <fct>             
 1 1     United States -10009207 Male   Age 0 to 4 years  
 2 1     United States -10322762 Male   Age 5 to 9 years  
 3 1     United States   9567476 Female Age 0 to 4 years  
 4 1     United States   7528626 Female Age 70 to 74 years
 5 1     United States   9873133 Female Age 5 to 9 years  
 6 1     United States -10618261 Male   Age 10 to 14 years
 7 1     United States  10180007 Female Age 10 to 14 years
 8 1     United States -11064752 Male   Age 20 to 24 years
 9 1     United States  10568188 Female Age 20 to 24 years
10 1     United States -12004570 Male   Age 25 to 29 years
# ℹ 26 more rows

Default:

usa_pyr |>
  ggplot(mapping = aes(x = value,
                       y = AGEGROUP,
                       fill = SEX)) +
  geom_col(width = 0.95, alpha = 0.75) +
  labs(x = "Count",
       y = NULL,
       title = "US Population Pyramid")

Fix the scales:

usa_pyr |>
  ggplot(mapping = aes(x = value,
                       y = AGEGROUP,
                       fill = SEX)) +
  geom_col(width = 0.95, alpha = 0.75) +
  # In the two scale functions we use \(x) label_number()
  # and \(x) str_remove_all() as anonymous functions.
  scale_x_continuous(
    labels = \(x) label_number(scale = 0.001, suffix = "k")(abs(x))) +
  scale_y_discrete(
    labels = \(x) str_remove_all(x, "Age | years")) +
  labs(x = "Count",
       y = NULL,
       title= "US Population Pyramid") +
  theme(legend.position = "bottom")

ARCOS Opioid Data

This example shows the results of setting up the ARCOS data as a duckdb database. Unlike most of the other examples we are not actually rendering this live on these pages because I do not want to include the >40GB ARCOS raw data file in the course repository.

First we take a look at a sample of the data; just the first thousand rows:

library(DBI) # will be autoloaded by dbplyr too
library(dbplyr)
library(here)

## Testing df
tmp <- read_tsv(here::here("raw", "arcos_all_washpost_noquotes.tsv"), n_max = 1000)
#> Rows: 1000 Columns: 16
#> ── Column specification ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────
#> Delimiter: "\t"
#> chr  (14): BUYER_DEA_NO, BUYER_BUS_ACT, BUYER_NAME, BUYER_ADDL_CO_INFO, BUYER_ADDRESS1, BUYER_ADDRESS2, BUYER_CITY, BUYER_COUNTY, REPOR...
#> dbl   (1): BUYER_ZIP
#> date  (1): date_iso
#> 
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

We’re going to import the date_iso column directly as a date by creating a custom class for it.

raw_colclasses <- rep("character", length(colnames(tmp)))
names(raw_colclasses) <- colnames(tmp)

## Hand-rolled date converter
setClass("myDate")
setAs("character", "myDate", function(from) as.Date(from, format="%Y-%m-%d"))


raw_colclasses[which(colnames(tmp) %in% "date_iso")] <- "myDate"
raw_colclasses
#>          BUYER_DEA_NO         BUYER_BUS_ACT            BUYER_NAME    BUYER_ADDL_CO_INFO        BUYER_ADDRESS1        BUYER_ADDRESS2 
#>           "character"           "character"           "character"           "character"           "character"           "character" 
#>            BUYER_CITY             BUYER_ZIP          BUYER_COUNTY         REPORTER_NAME              date_iso          Product_Name 
#>           "character"           "character"           "character"           "character"              "myDate"           "character" 
#> Combined_Labeler_Name  Revised_Company_Name       Reporter_family           BUYER_STATE 
#>  

When we import the data to duckdb in a moment we’ll use this raw_colclasses vector to set the column classes for duckdb’s duckdb_read_csv() function.

First we need to open a connection and say where the database file will exist.

# Raw data
path <- here::here("raw", "arcos_all_washpost_noquotes.tsv")

# DB connection
con <- dbConnect(duckdb::duckdb(),
                 dbdir = here("db", "arcos.duckdb"),
                 read_only = FALSE)

Now we read in the data using duckdb’s CSV reader. It’s parallelized and very fast. It will index all the columns by default.

## tictoc just to time things
tictoc::tic("Read DB")
duckdb::duckdb_read_csv(con, "arcos", path,
                        delim = "\t",
                        lower.case.names = TRUE,
                        colClasses = raw_colclasses)
#> Warning in duckdb::duckdb_read_csv(con, "arcos", path, delim = "\t", lower.case.names = TRUE, : Arguments passed to ... are currently not
#> used
tictoc::toc()
#> Read DB: 27.477 sec elapsed

This takes about 30 seconds on an MacBook Pro M3 Max with 128GB of RAM. On a Mac Mini with an M2 Pro processor and 16GB of RAM it takes about 80 seconds.

Now we can create a pseudo-table of our data. As a matter of convention we’ll use the _db prefix for objects that refer to the database but have not been materialized.


## Check
df_db <- tbl(con, "arcos")

df_db
#> # Source:   table<arcos> [?? x 16]
#> # Database: DuckDB v1.1.1 [root@Darwin 24.0.0:R 4.4.1//Users/kjhealy/Documents/data/arcos_opioids_test/db/arcos.duckdb]
#>    buyer_dea_no buyer_bus_act  buyer_name  buyer_addl_co_info buyer_address1 buyer_address2 buyer_city buyer_zip buyer_county reporter_name
#>    <chr>        <chr>          <chr>       <chr>              <chr>          <chr>          <chr>      <chr>     <chr>        <chr>        
#>  1 BC2511980    PRACTITIONER   COLEMAN, P… NA                 FAA C/O AMERI… PSC 801 BOX 41 LONDON     9498      NA           DIAMOND PHAR…
#>  2 BC2511980    PRACTITIONER   COLEMAN, P… NA                 FAA C/O AMERI… PSC 801 BOX 41 LONDON     9498      NA           DIAMOND PHAR…
#>  3 BC2462719    CHAIN PHARMACY CARR-GOTTS… CARRS PHARMACY #1… 1501 E HUFFMA… NA             ANCHORAGE  99511     ANCHORAGE    MCKESSON COR…
#>  4 AM3167120    PRACTITIONER   MERCHANT, … NA                 3710 WOODLAND… SUITE 1100     ANCHORAGE  99517     ANCHORAGE    DRX PHARMACE…
#>  5 AM3167120    PRACTITIONER   MERCHANT, … NA                 3710 WOODLAND… SUITE 1100     ANCHORAGE  99517     ANCHORAGE    DRX PHARMACE…
#>  6 BH6074431    PRACTITIONER   HALL, ROBE… NA                 3801 LAKE OTI… SUITE 300      ANCHORAGE  99508     ANCHORAGE    DRX PHARMACE…
#>  7 BH6074431    PRACTITIONER   HALL, ROBE… NA                 3801 LAKE OTI… SUITE 300      ANCHORAGE  99508     ANCHORAGE    DRX PHARMACE…
#>  8 BH6074431    PRACTITIONER   HALL, ROBE… NA                 3801 LAKE OTI… SUITE 300      ANCHORAGE  99508     ANCHORAGE    DRX PHARMACE…
#>  9 BJ5538371    PRACTITIONER   JOHNSTON, … ALASKA  SPINE INS… 3801 UNIVERSI… SUITE 300      ANCHORAGE  99508     ANCHORAGE    DRX PHARMACE…
#> 10 BM9540320    PRACTITIONER   MERCHANT, … NA                 1301 HUFFMAN … SUITE 205      ANCHORAGE  99515     ANCHORAGE    DRX PHARMACE…
#> # ℹ more rows
#> # ℹ 6 more variables: date_iso <date>, product_name <chr>, combined_labeler_name <chr>, revised_company_name <chr>, reporter_family <chr>,
#> #   buyer_state <chr>

Tabulating the Data

Now that our data are read in we can in many ways treat it just as if it existed as an ordinary tibble, even though it isn’t one.

How many rows do we have in our table anyway?

## How many rows?
df_db |> count() 
#> # Source:   SQL [1 x 1]
#> # Database: DuckDB v1.1.1 [root@Darwin 24.0.0:R 4.4.1//Users/kjhealy/Documents/data/arcos_opioids_test/db/arcos.duckdb]
#>           n
#>       <dbl>
#> 1 178598026

That’s 178,598,026 rows.

We can create new columns in the table. Not all functions are supported for this kind of operation.

## We could extract year and month from date using sql directly
# df_db <- df_db |>
#   mutate(year = sql("year(date_iso)"), 
#          month = sql("month(date_iso)"))  

## But we can now do the same thing using lubridate
## (Previously one couldn't do this)
df_db <- df_db |>
  mutate(year = lubridate::year(date_iso), 
         month = lubridate::month(date_iso))  
         
df_db
#> # Source:   SQL [?? x 18]
#> # Database: DuckDB v1.1.1 [root@Darwin 24.0.0:R 4.4.1//Users/kjhealy/Documents/data/arcos_opioids_test/db/arcos.duckdb]
#>    buyer_dea_no buyer_bus_act  buyer_name  buyer_addl_co_info buyer_address1 buyer_address2 buyer_city buyer_zip buyer_county reporter_name
#>    <chr>        <chr>          <chr>       <chr>              <chr>          <chr>          <chr>      <chr>     <chr>        <chr>        
#>  1 BC2511980    PRACTITIONER   COLEMAN, P… NA                 FAA C/O AMERI… PSC 801 BOX 41 LONDON     9498      NA           DIAMOND PHAR…
#>  2 BC2511980    PRACTITIONER   COLEMAN, P… NA                 FAA C/O AMERI… PSC 801 BOX 41 LONDON     9498      NA           DIAMOND PHAR…
#>  3 BC2462719    CHAIN PHARMACY CARR-GOTTS… CARRS PHARMACY #1… 1501 E HUFFMA… NA             ANCHORAGE  99511     ANCHORAGE    MCKESSON COR…
#>  4 AM3167120    PRACTITIONER   MERCHANT, … NA                 3710 WOODLAND… SUITE 1100     ANCHORAGE  99517     ANCHORAGE    DRX PHARMACE…
#>  5 AM3167120    PRACTITIONER   MERCHANT, … NA                 3710 WOODLAND… SUITE 1100     ANCHORAGE  99517     ANCHORAGE    DRX PHARMACE…
#>  6 BH6074431    PRACTITIONER   HALL, ROBE… NA                 3801 LAKE OTI… SUITE 300      ANCHORAGE  99508     ANCHORAGE    DRX PHARMACE…
#>  7 BH6074431    PRACTITIONER   HALL, ROBE… NA                 3801 LAKE OTI… SUITE 300      ANCHORAGE  99508     ANCHORAGE    DRX PHARMACE…
#>  8 BH6074431    PRACTITIONER   HALL, ROBE… NA                 3801 LAKE OTI… SUITE 300      ANCHORAGE  99508     ANCHORAGE    DRX PHARMACE…
#>  9 BJ5538371    PRACTITIONER   JOHNSTON, … ALASKA  SPINE INS… 3801 UNIVERSI… SUITE 300      ANCHORAGE  99508     ANCHORAGE    DRX PHARMACE…
#> 10 BM9540320    PRACTITIONER   MERCHANT, … NA                 1301 HUFFMAN … SUITE 205      ANCHORAGE  99515     ANCHORAGE    DRX PHARMACE…
#> # ℹ more rows
#> # ℹ 8 more variables: date_iso <date>, product_name <chr>, combined_labeler_name <chr>, revised_company_name <chr>, reporter_family <chr>,
#> #   buyer_state <chr>, year <dbl>, month <dbl>
         

Now we have year and month as their own columns.

We can try out counts and other summary operations without materializing them:

## Count, but not really, yet
tictoc::tic("Pretend counting")
by_state_yr_db <- df_db |>
  count(buyer_state, year, month)
tictoc::toc()
#> Pretend counting: 0.014 sec elapsed

## This object is still not real, 
## or "materialized" NB the ?? in the rows
## (the counts shown are correct though)
by_state_yr_db
#> # Source:   SQL [?? x 4]
#> # Database: DuckDB v1.1.1 [root@Darwin 24.0.0:R 4.4.1//Users/kjhealy/Documents/data/arcos_opioids_test/db/arcos.duckdb]
#>    buyer_state  year month     n
#>    <chr>       <dbl> <dbl> <dbl>
#>  1 AL           2012     7 58027
#>  2 AL           2007    12 40010
#>  3 AL           2010     6 53036
#>  4 AL           2012     4 54197
#>  5 AL           2010    10 47835
#>  6 AL           2007     1 39079
#>  7 AL           2008     5 42683
#>  8 AL           2011     5 52198
#>  9 AL           2006     7 36516
#> 10 AL           2006     4 34656
#> # ℹ more rows

And when we are ready we materialize the results with collect(). With data of this size, these operations are orders of magnitude faster than regular dplyr would be.

tictoc::tic("For reals counting")
by_state_yr <- by_state_yr_db |>
  collect()
tictoc::toc()
#> For reals counting: 0.175 sec elapsed

by_state_yr
#> # A tibble: 4,635 × 4
#>    buyer_state  year month     n
#>    <chr>       <dbl> <dbl> <dbl>
#>  1 AL           2010    12 53201
#>  2 AL           2009     3 44807
#>  3 AL           2009    11 43733
#>  4 AL           2011     4 49788
#>  5 AL           2010     1 41560
#>  6 AL           2012     5 59970
#>  7 AL           2011     7 49985
#>  8 AL           2008     4 42716
#>  9 AL           2008     7 45914
#> 10 AL           2007    10 43320
#> # ℹ 4,625 more rows

## Lots of data
by_state_yr |> 
  summarize(tot = sum(n))
#> # A tibble: 1 × 1
#>         tot
#>       <dbl>
#> 1 178598026