── 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
Example 09: Databases and APIs
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 ggplotLinking 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")- Get one from https://api.census.gov/data/key_signup.html
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 elapsedThis 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 178598026That’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 rowsAnd 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