Linking NaNDA With Your Data

Read-along notebook · October 7, 2026

Authors

Robert Melendez

You don’t need to run this

You don’t have to run this live to follow along. Robert Melendez runs the code on screen, and every output is already on this page. If you want to run it later, the code is in the materials ZIP.

Work in Stata or Python? Every code block has three tabs. The live session and the rendered output are R. The Stata and Python tabs show the same step in your language. Pick a tab once and the whole page follows.

What you’ll learn

By the end you will be able to:

Say which location identifier your data carry, how complete it is, and what it identifies. Section 1
Decide whether geocoding addresses is the right route for your study. Section 2
Tell a ZIP code from a ZCTA, and turn one into the other with a crosswalk without losing rows. Section 3
Test whether the rows you lost differ from the rows you kept, before anything from NaNDA is attached. Who did you lose?
Match the census boundary year across your codes, the crosswalk, and the NaNDA file, and recognize what goes wrong when they disagree. Section 4, side by side
Carry a measure past its last published year and write that decision into your methods. Carrying Social Services forward
Check a join: count who matched and name why the rest did not. Checking the join
Ask the four questions before every link, and cite the NaNDA datasets you use. Four questions, Section 5

0. Setup

library(tidyverse)   # dplyr for the joins, readr for the CSVs
library(readxl)      # the crosswalk is an .xlsx; installs with tidyverse, loads separately
library(here)        # file paths that work wherever this folder sits on your machine
library(janitor)     # tabyl(), for frequency tables you can read

* Run this from the workshop folder (the one holding workshop.Rproj), so the
* relative paths below resolve: File > Change Working Directory, or `cd`.
version 16
clear all
set more off

* One tempfile per table the R version keeps in memory. Stata holds one
* dataset at a time, so each step saves and the next step merges it back.
tempfile mydata socialservices parks ses2010 ses2020 xwalk mydata_zcta10 mydata_nanda

import pandas as pd            # data frames, read_csv / read_stata / merge
from pathlib import Path       # file paths that work wherever this folder sits

# The workshop folder: the one holding workshop.Rproj. Adjust if you moved the script.
ROOT = Path(__file__).resolve().parents[1] if "__file__" in globals() else Path.cwd()
DATA = ROOT / "data"

# Our own data. `zip` is a label: say so, or read_csv guesses "number" and
# 03042 becomes 3042.
mydata <- read_csv(
  here("data", "synthetic_data_v20260924.csv"),
  col_types = cols(zip = col_character())
)

# NaNDA: Social Services, ZCTA 2010, 1990-2022. Same rule for the ZCTA code.
socialservices_zcta10 <- read_csv(
  here("data", "nanda", "nanda_socials_Zcta10_1990-2022_01.csv"),
  col_types = cols(zcta10 = col_character())
)

# Parks and Socioeconomic Status come from ICPSR as .rda files. load() drops an
# object into your session under a name ICPSR picked, and returns that name.
# Print it, then rename it to something you will recognise later.
print(load(here("data", "nanda", "38586-0002-Data.rda")))

[1] "da38586.0002"

# The object keeps Robert's name, parks2022. The file is the 2018 ParkServe
# snapshot (DS0002, 2010 boundaries); ICPSR retitled it in October 2026.
parks2022_zcta10 <- da38586.0002
rm(da38586.0002)

print(load(here("data", "nanda", "38528-0003-Data.rda")))   # SES, ZCTA 2010, 2008-2017

[1] "da38528.0003"

ses2008_2017_zcta10 <- da38528.0003
rm(da38528.0003)

print(load(here("data", "nanda", "38528-0008-Data.rda")))   # SES, ZCTA 2020, 2018-2022

[1] "da38528.0008"

ses2018_2022_zcta20 <- da38528.0008
rm(da38528.0008)

* Our own data. Stata guesses "number" for a column of digits, and 03042 comes
* in as 3042. Check what it guessed and pad back to five characters if so.
import delimited using "data/synthetic_data_v20260924.csv", varnames(1) clear
capture confirm string variable zip
if _rc tostring zip, replace format(%05.0f)
* tostring writes a missing number as "."; blank it so missing() counts it.
replace zip = "" if zip == "."
save `mydata'

* NaNDA: Social Services, ZCTA 2010, 1990-2022. Same rule for the ZCTA code.
import delimited using "data/nanda/nanda_socials_Zcta10_1990-2022_01.csv", varnames(1) clear
capture confirm string variable zcta10
if _rc tostring zcta10, replace format(%05.0f)
save `socialservices'

* Parks and Socioeconomic Status: take the Stata format from ICPSR. Variable
* names arrive in upper case, as the ICPSR codebook lists them.
* The object keeps Robert's name, parks2022. The file is the 2018 ParkServe
* snapshot (DS0002, 2010 boundaries); ICPSR retitled it in October 2026.
use "data/nanda/38586-0002-Data.dta", clear          // Parks, ZCTA, 2018
capture confirm string variable ZCTA19
if _rc tostring ZCTA19, replace format(%05.0f)
save `parks'

use "data/nanda/38528-0003-Data.dta", clear          // SES, ZCTA 2010, 2008-2017
capture confirm string variable ZCTA10
if _rc tostring ZCTA10, replace format(%05.0f)
save `ses2010'

use "data/nanda/38528-0008-Data.dta", clear          // SES, ZCTA 2020, 2018-2022
capture confirm string variable ZCTA20
if _rc tostring ZCTA20, replace format(%05.0f)
save `ses2020'

# Our own data. `zip` is a label: say so, or pandas guesses "number" and
# 03042 becomes 3042.
mydata = pd.read_csv(DATA / "synthetic_data_v20260924.csv", dtype={"zip": str})

# NaNDA: Social Services, ZCTA 2010, 1990-2022. Same rule for the ZCTA code.
socialservices_zcta10 = pd.read_csv(
    DATA / "nanda" / "nanda_socials_Zcta10_1990-2022_01.csv", dtype={"zcta10": str}
)

# Parks and Socioeconomic Status: take the Stata format from ICPSR; pandas reads
# it without extra packages. Force every ZCTA code to five-character text.
def zcta_text(df, col):
    df[col] = df[col].astype(str).str.replace(r"\.0$", "", regex=True).str.zfill(5)
    return df

# The object keeps Robert's name, parks2022. The file is the 2018 ParkServe
# snapshot (DS0002, 2010 boundaries); ICPSR retitled it in October 2026.
parks2022_zcta10    = zcta_text(pd.read_stata(DATA / "nanda" / "38586-0002-Data.dta"), "ZCTA19")
ses2008_2017_zcta10 = zcta_text(pd.read_stata(DATA / "nanda" / "38528-0003-Data.dta"), "ZCTA10")
ses2018_2022_zcta20 = zcta_text(pd.read_stata(DATA / "nanda" / "38528-0008-Data.dta"), "ZCTA20")

Read ZIP codes as text, always

ZIP codes are labels. Read them as numbers and 03042 becomes 3042 with no warning, and every join for that row fails. That is all of New England, New Jersey, and Puerto Rico dropping out of your analysis.

glimpse(mydata)

Rows: 2,004
Columns: 11
$ sampleid            <dbl> 92202, 126800, 219095, 290998, 323367, 604662, 660…
$ year                <dbl> 2022, 2017, 2014, 2025, 2014, 2019, 2019, 2020, 20…
$ address             <chr> "5024 County Highway 33", NA, "MCDANIEL AVE", "118…
$ city                <chr> "Oneonta", "Eufaula", NA, "Brundidge", NA, NA, "Ma…
$ state               <chr> "AL", "AL", "AL", "AL", "AL", "AL", "AL", "AL", "A…
$ zip                 <chr> "35121", "36027", NA, "36010", NA, NA, "36555", NA…
$ social_interactions <dbl> 3, 1, 4, 1, 3, 5, 1, 1, 4, 4, 3, 3, 2, 3, 2, 2, 3,…
$ physical_activity   <dbl> 3, 3, 2, 4, 3, 3, 2, 2, 1, 3, 2, 3, 4, 2, 4, 2, 3,…
$ hot_meal_days       <dbl> 7, 7, 7, 6, 6, 5, 4, 6, 7, 7, 6, 3, 5, 5, 5, 4, 7,…
$ loneliness          <dbl> 1.3150288, 0.9892247, 0.1722808, 1.4024957, 2.9253…
$ own_pet             <dbl> 1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 1, 0,…

use `mydata', clear
describe
list in 1/5

mydata.info()
mydata.head()

1. Which identifiers your data carry

Most survey and clinical datasets do not have geography. They have a string: an address, a city, a ZIP code, collected for mailing or billing, not for placing people in neighborhoods. That is the job now: to place each person in a neighborhood that another dataset has measured.

Before any merge, answer three questions about that string:

Which identifier do I have? An address, a ZIP code, a Census code?
How complete is it? Rows with nothing to look up never match.
What is it an identifier of? A ZIP code is a set of mail routes. An address becomes an identifier only after geocoding.

If your data already carry Census tract or ZCTA codes from an earlier geocoding step, those codes follow the boundaries of one census year, and the boundary question in Section 4 applies to them too.

What’s in the synthetic dataset

Variable	Type	Notes
`sampleid`	numeric	2,004 unique; one row per person
`year`	numeric	2014–2025
`address`, `city`, `state`, `zip`	character	`zip` missing for ~17% of rows
`social_interactions`, `physical_activity`	1–5
`hot_meal_days`	1–7
`loneliness`	continuous
`own_pet`	0/1

What NaNDA is published at

NaNDA describes neighborhoods, not people. Every measure is a count, a density, or a rate for a geographic unit in a given year. There are no individual records in the archive, so the linking is on you. NaNDA has no record of where your participants live.

Today we use three datasets from the NaNDA catalog on ICPSR:

NaNDA: Social Services (1990–2022)
NaNDA: Parks (2018, 2022)
NaNDA: Socioeconomic Status (2008–2022, two files)

All three are released at the Census tract level and the ZCTA level. We use the ZCTA files, because ZIP codes are what the synthetic data can offer. Every file we join to is drawn on 2010 boundaries. Robert Melendez explains that choice in the session, and Section 4 shows what changes when you mix boundary years.

Our file has ZIP codes and NaNDA has ZCTAs. They look alike, usually the same five digits, and they are not the same thing. Section 3 is about that gap.

# Count the rows that have a ZIP.
mydata %>%
  summarise(
    rows        = n(),
    has_zip     = sum(!is.na(zip)),
    missing_zip = sum(is.na(zip)),
    pct_missing = round(100 * mean(is.na(zip)), 1)
  )

rows	has_zip	missing_zip	pct_missing
2004	1662	342	17.1

* Count the rows that have a ZIP.
use `mydata', clear
generate byte missing_zip = missing(zip)
summarize missing_zip
display "rows = " r(N) ", missing ZIP = " r(sum) ", pct missing = " %4.1f 100 * r(mean)

# Count the rows that have a ZIP.
pd.Series({
    "rows":        len(mydata),
    "has_zip":     mydata["zip"].notna().sum(),
    "missing_zip": mydata["zip"].isna().sum(),
    "pct_missing": round(100 * mydata["zip"].isna().mean(), 1),
})

# address, city and zip are missing on different rows. A row with no ZIP may
# still have an address you could geocode, and the other way round, so the
# usable geography depends on the route.
mydata %>%
  count(
    address_missing = is.na(address),
    city_missing    = is.na(city),
    zip_missing     = is.na(zip)
  ) %>%
  arrange(desc(n))

address_missing	city_missing	zip_missing	n
FALSE	FALSE	FALSE	1618
TRUE	TRUE	TRUE	220
FALSE	TRUE	TRUE	105
TRUE	FALSE	FALSE	26
FALSE	TRUE	FALSE	15
FALSE	FALSE	TRUE	13
TRUE	FALSE	TRUE	4
TRUE	TRUE	FALSE	3

* address, city and zip are missing on different rows. A row with no ZIP may
* still have an address you could geocode, and the other way round, so the
* usable geography depends on the route.
use `mydata', clear
generate byte address_missing = missing(address)
generate byte city_missing    = missing(city)
generate byte zip_missing     = missing(zip)
contract address_missing city_missing zip_missing, freq(n)
gsort -n
list, clean noobs

# address, city and zip are missing on different rows. A row with no ZIP may
# still have an address you could geocode, and the other way round, so the
# usable geography depends on the route.
(mydata
   .assign(address_missing=mydata["address"].isna(),
           city_missing=mydata["city"].isna(),
           zip_missing=mydata["zip"].isna())
   .value_counts(["address_missing", "city_missing", "zip_missing"])
   .rename("n")
   .reset_index()
   .sort_values("n", ascending=False))

2. If you have addresses: geocoding

A street address is the most precise location most studies collect, but it cannot be joined to anything until it is geocoded. Geocoding turns it into a point on the map. A spatial join then places that point inside a Census tract, and the tract-level NaNDA file joins on the tract code.

The route, step by step:

Geocode. Send each address to a geocoding service and get back coordinates plus a precision code. The Census Bureau’s geocoder is free and returns the tract directly. tidygeocoder in R talks to OpenStreetMap, ArcGIS, Google, and others through one interface.
Check precision. A rooftop match and a “somewhere in this city” match both come back with coordinates. Keep the tract only where the match is precise enough to trust. With OpenStreetMap, that is the place rank.
Spatial join to 2010 tracts. Download the TIGER/Line tract boundaries for the census year your NaNDA file uses and intersect the points with them.
Join to the tract-level NaNDA file on the tract code, then run the same diagnostics as Section 4.

Why we are not doing it live:

It is slow and it costs. Free services throttle to about one address per second. Paid ones charge per thousand addresses.
It is a terms-of-service question. Some services restrict storing or sharing the results, so a geocoded file may not be shareable even when the addresses are.
It fails differently. Addresses that will not parse, apartment-level ambiguity, and PO boxes all produce a coordinate or nothing, and the failure never looks like a join error.

The shape of it in R, using OpenStreetMap:

library(tidygeocoder)   # one interface to many geocoding services
library(tigris)         # Census boundary files
library(sf)             # spatial joins

mydata_geocoded <- mydata %>%
  geocode(street = address, city = city, state = state, postalcode = zip,
          method = "osm", full_results = TRUE)

# Keep a tract only where the match is precise enough (OSM place_rank >= 26).
points <- mydata_geocoded %>%
  filter(!is.na(lat)) %>%
  st_as_sf(coords = c("long", "lat"), crs = 4326)

tracts2010 <- tracts(state = unique(mydata$state), cb = FALSE, year = 2010)
points <- st_transform(points, st_crs(tracts2010))

mydata_tract10 <- st_join(points, tracts2010, join = st_intersects) %>%
  mutate(tract_fips10 = if_else(place_rank >= 26, GEOID10, NA_character_)) %>%
  st_drop_geometry()

From here the tract-level NaNDA files join on tract_fips10 exactly as the ZCTA files join on zcta10 below. Everything in Section 4 applies unchanged. The complete script is in your materials as r/geocode_to_tract.R, with a Python version beside it and a Stata version that goes through the Census Geocoder instead; each language folder’s README says what it needs and how long it takes.

3. ZIP is not ZCTA

A ZIP code is a Postal Service artifact. It names a set of delivery routes, and it was built for mail delivery.

Some ZIP codes are a single building.
Some are a bank of post office boxes with no residents.
Some cross county and state lines.
The Postal Service redraws them whenever the mail requires, without telling anyone doing research.

A ZCTA, a ZIP Code Tabulation Area, is the Census Bureau’s answer to the fact that researchers use ZIP codes anyway. For each decennial census, the Bureau assigns every census block to the ZIP code most of its addresses carry. The result is a set of areas with the same five-digit names as ZIP codes that tile the country and can be measured. NaNDA’s ZCTA files are built on these.

Most of the time a ZIP code and the ZCTA of the same name cover roughly the same people. ZIP codes with no ZCTA, such as post office boxes and single-site ZIPs, have nowhere to go. ZIPs mostly absorbed into a neighboring ZCTA are assigned there, under a different number. So the step from ZIP to ZCTA is a lookup against a crosswalk, not a rename of the column. And the crosswalk follows the boundaries of one census year, like everything else.

The crosswalk we’re using

UDS Mapper ZIP Code to ZCTA Crosswalk (2019 release, archived copy). Included in the ZIP as data/zip_to_zcta_2019.xlsx.
NaNDA’s published merge code: Chenoweth & Khan, code for merging NaNDA ZCTA-level datasets with the UDS Mapper crosswalk (2024).

Use the published merge code rather than a snippet copied from a screen. It is citable and maintained.

# Every column as text, for the same reason as `zip` above. Then pad back to five
# characters: whatever stripped a leading zero upstream did it before we got here.
zipzctaxwalk2019 <- read_excel(
  here("data", "zip_to_zcta_2019.xlsx"),
  col_types = "text"
) %>%
  mutate(
    ZIP_CODE = str_pad(ZIP_CODE, 5, pad = "0"),
    ZCTA     = str_pad(ZCTA,     5, pad = "0")
  )

# zip_join_type records how each ZIP was assigned its ZCTA, or that it has
# none. Check it before joining.
tabyl(zipzctaxwalk2019, zip_join_type)

zip_join_type	n	percent
Spatial join to ZCTA	7960	0.1936410
Territory ZCTA Add	2	0.0000487
Zip Matches ZCTA	33120	0.8057022
populated ZCTA, missing zip	16	0.0003892
territory zip, no ZCTA available	9	0.0002189

* Every column as text, for the same reason as `zip` above. Then pad back to
* five characters: whatever stripped a leading zero upstream did it before we
* got here.
import excel using "data/zip_to_zcta_2019.xlsx", firstrow allstring clear
replace ZIP_CODE = substr("00000", 1, 5 - strlen(ZIP_CODE)) + ZIP_CODE if strlen(ZIP_CODE) < 5
replace ZCTA     = substr("00000", 1, 5 - strlen(ZCTA))     + ZCTA     if strlen(ZCTA)     < 5
keep ZIP_CODE ZCTA zip_join_type
save `xwalk'

* zip_join_type records how each ZIP was assigned its ZCTA, or that it has
* none. Check it before joining.
tabulate zip_join_type

# Every column as text, for the same reason as `zip` above. Then pad back to
# five characters: whatever stripped a leading zero upstream did it before we
# got here. (read_excel needs the openpyxl package: pip install openpyxl)
zipzctaxwalk2019 = pd.read_excel(DATA / "zip_to_zcta_2019.xlsx", dtype=str)
zipzctaxwalk2019["ZIP_CODE"] = zipzctaxwalk2019["ZIP_CODE"].str.zfill(5)
zipzctaxwalk2019["ZCTA"]     = zipzctaxwalk2019["ZCTA"].str.zfill(5)

# zip_join_type records how each ZIP was assigned its ZCTA, or that it has
# none. Check it before joining.
zipzctaxwalk2019["zip_join_type"].value_counts(dropna=False)

mydata_zcta10 <- mydata %>%
  left_join(
    zipzctaxwalk2019 %>% select(ZIP_CODE, ZCTA, zip_join_type),
    by = c("zip" = "ZIP_CODE")
  ) %>%
  rename(zcta10 = ZCTA, zip_join_type10 = zip_join_type)

# A left join keeps every row of the left-hand table. Compare the counts.
# If this number grew, the crosswalk has more than one row per ZIP and some
# people are now counted twice.
c(rows_before = nrow(mydata), rows_after = nrow(mydata_zcta10))

rows_before  rows_after 
       2004        2004

# A row can lack a ZCTA for two reasons. A row with no ZIP has nothing to look
# up. A row whose ZIP is missing from the crosswalk had a ZIP but no match: PO-box ZIPs, single-building ZIPs, ZIPs the USPS has retired.
mydata_zcta10 %>%
  mutate(zcta_status = case_when(
    is.na(zip)    ~ "no ZIP",
    is.na(zcta10) ~ "ZIP not in crosswalk",
    TRUE          ~ "has a ZCTA"
  )) %>%
  tabyl(zcta_status)

zcta_status	n	percent
ZIP not in crosswalk	2	0.0009980
has a ZCTA	1660	0.8283433
no ZIP	342	0.1706587

use `xwalk', clear
rename (ZIP_CODE ZCTA zip_join_type) (zip zcta10 zip_join_type10)
save `xwalk', replace

use `mydata', clear
count
local rows_before = r(N)

* m:1 means many rows of ours to one row of theirs. If the crosswalk has two
* rows for one ZIP, the merge stops with an error. keep(master match) is the
* left join: every row of ours stays, matched or not.
merge m:1 zip using `xwalk', keep(master match)
count
display "rows before = `rows_before', rows after = " r(N)

* A row can lack a ZCTA for two reasons. A row with no ZIP has nothing to look
* up. A row whose ZIP is missing from the crosswalk had a ZIP but no match: PO-box ZIPs, single-building ZIPs, retired ZIPs.
generate str zcta_status = "has a ZCTA"
replace  zcta_status = "ZIP not in crosswalk" if _merge == 1
replace  zcta_status = "no ZIP"               if missing(zip)
tabulate zcta_status
drop _merge
save `mydata_zcta10'

mydata_zcta10 = (
    mydata.merge(
        zipzctaxwalk2019[["ZIP_CODE", "ZCTA", "zip_join_type"]],
        how="left", left_on="zip", right_on="ZIP_CODE",
    )
    .drop(columns="ZIP_CODE")
    .rename(columns={"ZCTA": "zcta10", "zip_join_type": "zip_join_type10"})
)

# A left join keeps every row of the left-hand table. Compare the counts.
# If this number grew, the crosswalk has more than one row per ZIP and some
# people are now counted twice.
{"rows_before": len(mydata), "rows_after": len(mydata_zcta10)}

# A row can lack a ZCTA for two reasons. A row with no ZIP has nothing to look
# up. A row whose ZIP is missing from the crosswalk had a ZIP but no match: PO-box ZIPs, single-building ZIPs, retired ZIPs.
import numpy as np
mydata_zcta10.assign(zcta_status=np.select(
    [mydata_zcta10["zip"].isna(), mydata_zcta10["zcta10"].isna()],
    ["no ZIP", "ZIP not in crosswalk"], default="has a ZCTA",
))["zcta_status"].value_counts()

Which ZCTA boundaries does your crosswalk produce?

The 2019 release of the UDS Mapper ZIP Code to ZCTA Crosswalk produces 2010-boundary ZCTAs. That makes it the right partner for NaNDA’s 2010-boundary files, which is what we use throughout. With your own data, read the crosswalk’s documentation for its target boundaries before anything else.

Who did you lose?

The table above gives the size of the loss. The diagnostic is whether the rows that did not get a ZCTA differ from the rows that did, and it belongs here, before anything from NaNDA is attached.

In this file about 17% of rows have no ZIP at all. If those rows are disproportionately rural, older, or in worse health, the matched sample is no longer the sample you drew. Every neighborhood coefficient you estimate is then conditional on having a mailable address. That is a selection problem, and no output flags it.

Neither language does this for you. Isolate the rows without a ZCTA and put their means beside the others’ means on the variables the study is about.

# Compare the rows that matched with those that did not, on the variables the
# study is about.
mydata_zcta10 %>%
  group_by(has_zcta = !is.na(zcta10)) %>%
  summarise(
    n                 = n(),
    loneliness        = round(mean(loneliness,        na.rm = TRUE), 2),
    physical_activity = round(mean(physical_activity, na.rm = TRUE), 2),
    own_pet           = round(mean(own_pet,           na.rm = TRUE), 3),
    .groups = "drop"
  )

has_zcta	n	loneliness	physical_activity	own_pet
FALSE	344	1.14	3.00	0.660
TRUE	1660	1.11	3.04	0.638

* Compare the rows that matched with those that did not, on the variables the
* study is about.
use `mydata_zcta10', clear
generate byte has_zcta = !missing(zcta10)
tabstat loneliness physical_activity own_pet, by(has_zcta) statistics(n mean) format(%9.3f)

# Compare the rows that matched with those that did not, on the variables the
# study is about.
(mydata_zcta10.assign(has_zcta=mydata_zcta10["zcta10"].notna())
   .groupby("has_zcta")
   .agg(n=("has_zcta", "size"),
        loneliness=("loneliness", "mean"),
        physical_activity=("physical_activity", "mean"),
        own_pet=("own_pet", "mean"))
   .round({"loneliness": 2, "physical_activity": 2, "own_pet": 3})
   .reset_index())

The same check applies after every join that follows. Section 4 repeats it against the NaNDA file, where a second set of reasons can drop rows.

4. Which boundaries? Joining, and checking the join

The join itself rarely goes wrong. If a column name is misspelled the software says so. The errors come from the inputs and from what you do with the result, and they are easy to miss.

Silent error 1: boundaries. ZCTA boundaries were redrawn for 2020. Most codes kept their five digits but changed their edges; some were added, dropped, or split. A 2010 code joined to a 2020 file usually finds a row, and the row usually describes approximately the same place. Your results will not show which matches were approximate.
Silent error 2: never looking at what fell out. A left join keeps every row of your data and fills NaNDA’s columns with NA where nothing matched. An inner join drops those rows. Either way the analysis proceeds. Stata’s merge reports _merge. R’s dplyr joins report nothing, so the count has to be built by hand.

The rest of this section is mostly diagnostics.

Carrying Social Services forward

Social Services ends in 2022; our data runs to 2025. Before joining, we give every ZCTA rows through 2025 and carry its most recent values into them. This is a modeling decision: the 2023 to 2025 rows will match, but they carry 2022’s numbers. State it in your methods.

# complete() adds the missing ZCTA-year rows; fill() carries each ZCTA's last
# value down into them. fill() does not stop at 2022: any value missing inside
# 1990-2022 (population, and the densities where population is zero) also gets
# the year before.
socialservices_zcta10 <- socialservices_zcta10 %>%
  complete(zcta10, year = 1990:2025) %>%
  arrange(zcta10, year) %>%
  group_by(zcta10) %>%
  fill(everything(), .direction = "down") %>%
  ungroup()

# Every year from 2022 on should now carry the same values.
tabyl(socialservices_zcta10, year) %>% filter(year >= 2020)

year	n	percent
2020	33144	0.0277778
2021	33144	0.0277778
2022	33144	0.0277778
2023	33144	0.0277778
2024	33144	0.0277778
2025	33144	0.0277778

* Stata has no complete(); every ZCTA here has a 2022 row, so make three copies
* of it and number them 2023-2025. Then carry each ZCTA's last value down into
* anything still missing. Like R's fill(), that also reaches values missing
* inside 1990-2022 (population, and densities where population is zero).
use `socialservices', clear
expand 4 if year == 2022
bysort zcta10 year: replace year = 2021 + _n if year == 2022
sort zcta10 year
foreach v of varlist _all {
    if !inlist("`v'", "zcta10", "year") {
        by zcta10: replace `v' = `v'[_n-1] if missing(`v')
    }
}
save `socialservices', replace

* Every year from 2022 on should now carry the same values.
tabulate year if year >= 2020

# pandas has no complete(); every ZCTA here has a 2022 row, so copy it to
# 2023-2025. Then carry each ZCTA's last value down into anything still missing.
# Like R's fill(), that also reaches values missing inside 1990-2022
# (population, and densities where population is zero).
ss_2022 = socialservices_zcta10[socialservices_zcta10["year"] == 2022]
socialservices_zcta10 = (
    pd.concat([socialservices_zcta10] + [ss_2022.assign(year=y) for y in (2023, 2024, 2025)])
    .sort_values(["zcta10", "year"])
    .reset_index(drop=True)
)
cols = socialservices_zcta10.columns.drop("zcta10")
socialservices_zcta10[cols] = socialservices_zcta10.groupby("zcta10")[cols].ffill()

# Every year from 2022 on should now carry the same values.
socialservices_zcta10.loc[socialservices_zcta10["year"] >= 2020, "year"].value_counts().sort_index()

The join

mydata_nanda <- mydata_zcta10 %>%
  # Social Services is longitudinal, so the key is ZCTA *and* year.
  # `in_nanda` is added before the join: it is TRUE on every row of the NaNDA
  # side, so after the join it is TRUE where something matched and NA where
  # nothing did. Stata's `_merge` flags matched rows; here `in_nanda` does the
  # same.
  left_join(
    socialservices_zcta10 %>%
      select(zcta10, year, totpop, aland10,
             count_totindivfamilyservices, den_totindivfamilyservices,
             count_childyouthservices,     den_childyouthservices,
             count_elderservices,          den_elderservices) %>%
      mutate(in_nanda = TRUE),
    by = c("zcta10", "year")
  ) %>%
  mutate(in_nanda = coalesce(in_nanda, FALSE)) %>%
  # Parks is a single 2018 snapshot, so it joins on ZCTA alone and every year of
  # a person's records gets the same value. That assumes park provision held
  # still across the study period.
  left_join(
    parks2022_zcta10 %>%
      select(ZCTA19, ANY_OPEN_PARK, COUNT_OPEN_PARKS_TC10, PROP_PARK_AREA_ZCTA),
    by = c("zcta10" = "ZCTA19")
  )

c(rows_before = nrow(mydata_zcta10), rows_after = nrow(mydata_nanda))

rows_before  rows_after 
       2004        2004

use `socialservices', clear
keep zcta10 year totpop aland10 ///
     count_totindivfamilyservices den_totindivfamilyservices ///
     count_childyouthservices     den_childyouthservices ///
     count_elderservices          den_elderservices
save `socialservices', replace

use `mydata_zcta10', clear
count
local rows_before = r(N)

* Social Services is longitudinal, so the key is ZCTA *and* year. Stata's
* _merge is the diagnostic the R version has to build by hand: 3 = matched,
* 1 = ours only. Keep it as in_nanda.
merge m:1 zcta10 year using `socialservices', keep(master match)
generate byte in_nanda = (_merge == 3)
drop _merge

* Parks is a single 2018 snapshot, so it joins on ZCTA alone and every year of
* a person's records gets the same value. That assumes park provision held
* still across the study period.
preserve
    use `parks', clear
    keep ZCTA19 ANY_OPEN_PARK COUNT_OPEN_PARKS_TC10 PROP_PARK_AREA_ZCTA
    rename ZCTA19 zcta10
    save `parks', replace
restore
merge m:1 zcta10 using `parks', keep(master match) nogenerate

count
display "rows before = `rows_before', rows after = " r(N)
save `mydata_nanda'

# Social Services is longitudinal, so the key is ZCTA *and* year.
# merge(indicator=True) does what Stata's _merge does: "both"
# where something matched, "left_only" where nothing did. Keep it as in_nanda.
mydata_nanda = (
    mydata_zcta10.merge(
        socialservices_zcta10[[
            "zcta10", "year", "totpop", "aland10",
            "count_totindivfamilyservices", "den_totindivfamilyservices",
            "count_childyouthservices",     "den_childyouthservices",
            "count_elderservices",          "den_elderservices",
        ]],
        how="left", on=["zcta10", "year"], indicator=True,
    )
    .assign(in_nanda=lambda d: d["_merge"].eq("both"))
    .drop(columns="_merge")
    # Parks is a single 2018 snapshot, so it joins on ZCTA alone and every year
    # of a person's records gets the same value. That assumes park provision
    # held still across the study period.
    .merge(
        parks2022_zcta10[["ZCTA19", "ANY_OPEN_PARK", "COUNT_OPEN_PARKS_TC10", "PROP_PARK_AREA_ZCTA"]],
        how="left", left_on="zcta10", right_on="ZCTA19",
    )
    .drop(columns="ZCTA19")
)

{"rows_before": len(mydata_zcta10), "rows_after": len(mydata_nanda)}

Socioeconomic status, on 2010 boundaries

NaNDA’s SES measures are five-year ACS estimates, so each file’s value stands for a stretch of years. NaNDA publishes two ZCTA files: 2008–2017 on 2010 boundaries and 2018–2022 on 2020 boundaries. Our crosswalk produces 2010-boundary codes, so every row joins the 2010-boundary file. The 2020-boundary file appears once, below, to show what changes if you mix them.

Limitation: affluence for observations from 2018 onward is the 2013–2017 estimate, because NaNDA publishes no 2010-boundary ZCTA file for later years.

# The 2010-boundary file (DS0003) names its affluence measure for the ACS years
# behind it. Give it a plain name before joining.
mydata_nanda <- mydata_nanda %>%
  left_join(
    ses2008_2017_zcta10 %>%
      select(ZCTA10, AFFLUENCE13_17) %>%
      rename(AFFLUENCE = AFFLUENCE13_17),
    by = c("zcta10" = "ZCTA10")
  )

c(rows = nrow(mydata_nanda), with_affluence = sum(!is.na(mydata_nanda$AFFLUENCE)))

          rows with_affluence 
          2004           1657

* The 2010-boundary file (DS0003) names its affluence measure for the ACS years
* behind it. Give it a plain name before merging. The 2020-boundary file is
* trimmed the same way here, for the comparison below only.
use `ses2010', clear
keep ZCTA10 AFFLUENCE13_17
rename (ZCTA10 AFFLUENCE13_17) (zcta10 AFFLUENCE)
save `ses2010', replace

use `ses2020', clear
keep ZCTA20 AFFLUENCE
rename ZCTA20 zcta10
save `ses2020', replace

use `mydata_nanda', clear
merge m:1 zcta10 using `ses2010', keep(master match) nogenerate
count if !missing(AFFLUENCE)
display "rows = " _N ", with affluence = " r(N)
save `mydata_nanda', replace

# The 2010-boundary file (DS0003) names its affluence measure for the ACS years
# behind it. Give it a plain name before merging. The 2020-boundary file is
# trimmed the same way here, for the comparison below only.
ses2010 = ses2008_2017_zcta10[["ZCTA10", "AFFLUENCE13_17"]].rename(
    columns={"ZCTA10": "zcta10", "AFFLUENCE13_17": "AFFLUENCE"})
ses2020 = ses2018_2022_zcta20[["ZCTA20", "AFFLUENCE"]].rename(columns={"ZCTA20": "zcta10"})

mydata_nanda = mydata_nanda.merge(ses2010, how="left", on="zcta10")

{"rows": len(mydata_nanda), "with_affluence": mydata_nanda["AFFLUENCE"].notna().sum()}

2010 vs 2020 boundaries, side by side

The same join, run twice with the same crosswalked codes: once against the 2010-boundary SES file and once against the 2020-boundary file. Both runs match the same rows at the same rate, because most ZCTAs kept their five digits from one census to the next. The value attached to each row differs: the same people get a different mean affluence from each file. The output does not show which run used matching boundaries; that comes from the documentation of each input.

# Socioeconomic Status, 2010-boundary file (ICPSR 38528, DS0003), joined to
# every row. Same rename as above: the two files call the same construct
# different things.
on_2010 <- mydata_zcta10 %>%
  left_join(
    ses2008_2017_zcta10 %>%
      select(ZCTA10, AFFLUENCE13_17) %>%
      rename(AFFLUENCE = AFFLUENCE13_17),
    by = c("zcta10" = "ZCTA10")
  )

on_2010 %>%
  summarise(
    rows           = n(),
    matched        = sum(!is.na(AFFLUENCE)),
    match_rate_pct = round(100 * mean(!is.na(AFFLUENCE)), 1),
    mean_affluence = round(mean(AFFLUENCE, na.rm = TRUE), 3)
  )

rows	matched	match_rate_pct	mean_affluence
2004	1657	82.7	0.382

* Socioeconomic Status, 2010-boundary file (ICPSR 38528, DS0003), already
* trimmed and renamed above, against every row.
use `mydata_zcta10', clear
merge m:1 zcta10 using `ses2010', keep(master match) nogenerate
generate byte matched = !missing(AFFLUENCE)
summarize matched
local rows = r(N)
local matched = r(sum)
local rate = 100 * r(mean)
summarize AFFLUENCE
display "SES, ZCTA 2010 file: rows = `rows', matched = `matched', match rate = " ///
        %4.1f `rate' "%, mean affluence = " %6.3f r(mean)
tempfile on_2010
save `on_2010'

# Socioeconomic Status, 2010-boundary file (ICPSR 38528, DS0003), joined to
# every row. Same rename as above: the two files call the same construct
# different things.
on_2010 = mydata_zcta10.merge(
    ses2008_2017_zcta10[["ZCTA10", "AFFLUENCE13_17"]].rename(columns={"AFFLUENCE13_17": "AFFLUENCE"}),
    how="left", left_on="zcta10", right_on="ZCTA10",
).drop(columns="ZCTA10")

def match_summary(df):
    return pd.Series({
        "rows":           len(df),
        "matched":        df["AFFLUENCE"].notna().sum(),
        "match_rate_pct": round(100 * df["AFFLUENCE"].notna().mean(), 1),
        "mean_affluence": round(df["AFFLUENCE"].mean(), 3),
    })

match_summary(on_2010)

# The same 2010-boundary codes, straight from the crosswalk, against the
# 2020-boundary file (ICPSR 38528, DS0008). The codes and the boundaries come
# from different censuses, and the join still returns a match rate and a mean.
on_2020 <- mydata_zcta10 %>%
  left_join(
    ses2018_2022_zcta20 %>% select(ZCTA20, AFFLUENCE),
    by = c("zcta10" = "ZCTA20")
  )

on_2020 %>%
  summarise(
    rows           = n(),
    matched        = sum(!is.na(AFFLUENCE)),
    match_rate_pct = round(100 * mean(!is.na(AFFLUENCE)), 1),
    mean_affluence = round(mean(AFFLUENCE, na.rm = TRUE), 3)
  )

rows	matched	match_rate_pct	mean_affluence
2004	1657	82.7	0.36

* The same 2010-boundary codes, straight from the crosswalk, against the
* 2020-boundary file (ICPSR 38528, DS0008). The codes and the boundaries come
* from different censuses, and the merge still returns a match rate and a mean.
use `mydata_zcta10', clear
merge m:1 zcta10 using `ses2020', keep(master match) nogenerate
generate byte matched = !missing(AFFLUENCE)
summarize matched
local rows = r(N)
local matched = r(sum)
local rate = 100 * r(mean)
summarize AFFLUENCE
display "SES, ZCTA 2020 file: rows = `rows', matched = `matched', match rate = " ///
        %4.1f `rate' "%, mean affluence = " %6.3f r(mean)
tempfile on_2020
save `on_2020'

# The same 2010-boundary codes, straight from the crosswalk, against the
# 2020-boundary file (ICPSR 38528, DS0008). The codes and the boundaries come
# from different censuses, and the join still returns a match rate and a mean.
on_2020 = mydata_zcta10.merge(
    ses2018_2022_zcta20[["ZCTA20", "AFFLUENCE"]],
    how="left", left_on="zcta10", right_on="ZCTA20",
).drop(columns="ZCTA20")

match_summary(on_2020)

bind_rows(
  `SES, ZCTA 2010 file` = on_2010,
  `SES, ZCTA 2020 file` = on_2020,
  .id = "nanda_file"
) %>%
  group_by(nanda_file) %>%
  summarise(
    rows           = n(),
    matched        = sum(!is.na(AFFLUENCE)),
    match_rate_pct = round(100 * mean(!is.na(AFFLUENCE)), 1),
    mean_affluence = round(mean(AFFLUENCE, na.rm = TRUE), 3),
    .groups = "drop"
  )

nanda_file	rows	matched	match_rate_pct	mean_affluence
SES, ZCTA 2010 file	2004	1657	82.7	0.382
SES, ZCTA 2020 file	2004	1657	82.7	0.360

use `on_2010', clear
generate str nanda_file = "SES, ZCTA 2010 file"
append using `on_2020'
replace nanda_file = "SES, ZCTA 2020 file" if missing(nanda_file)
generate byte one = 1
collapse (sum) rows = one (sum) matched (mean) match_rate_pct = matched ///
         (mean) mean_affluence = AFFLUENCE, by(nanda_file)
replace match_rate_pct = round(100 * match_rate_pct, 0.1)
replace mean_affluence = round(mean_affluence, 0.001)
list, clean noobs

pd.concat(
    {"SES, ZCTA 2010 file": on_2010, "SES, ZCTA 2020 file": on_2020},
    names=["nanda_file"],
).groupby(level="nanda_file").apply(match_summary)

The match rate is identical in both rows, and the mean differs. A boundary error shows up as a wrong value on a matched row, so a match rate cannot detect it.

Match the boundaries, and say which ones you matched

Every code in a join follows the boundaries of one census year: the codes in your file, the codes the crosswalk produces, and the codes the NaNDA file is drawn on. All three have to agree. Ours agree on 2010. If your data run past 2019 and you want 2020 boundaries instead, use a crosswalk that produces 2020 codes and the 2020-boundary NaNDA files, and state the choice in your methods. Choose the boundary year before the join.

Checking the join

The same question as Section 3, one join later. A row that reached this point has a ZCTA; it can still fail to find a NaNDA row. dplyr’s joins report nothing, and Stata’s _merge gives counts but not reasons, so we build both: the match rate, then the reason for every miss.

# Overall match rate.
mydata_nanda %>%
  summarise(
    rows           = n(),
    matched        = sum(in_nanda),
    match_rate_pct = round(100 * mean(in_nanda), 1)
  )

rows	matched	match_rate_pct
2004	1658	82.7

# The same number, by year. 2023 onward matches because we carried Social
# Services forward; without that step these rows would show 0%.
mydata_nanda %>%
  group_by(year) %>%
  summarise(
    rows           = n(),
    matched        = sum(in_nanda),
    match_rate_pct = round(100 * mean(in_nanda), 1),
    .groups = "drop"
  )

year	rows	matched	match_rate_pct
2014	167	132	79.0
2015	182	146	80.2
2016	160	127	79.4
2017	200	161	80.5
2018	164	135	82.3
2019	164	139	84.8
2020	157	137	87.3
2021	168	136	81.0
2022	167	138	82.6
2023	161	137	85.1
2024	154	132	85.7
2025	160	138	86.2

* Overall match rate.
use `mydata_nanda', clear
summarize in_nanda
display "rows = " r(N) ", matched = " r(sum) ", match rate = " %4.1f 100 * r(mean) "%"

* The same number, by year. 2023 onward matches because we carried Social
* Services forward; without that step these rows would show 0%.
tabulate year in_nanda, row nofreq

# Overall match rate.
pd.Series({
    "rows":           len(mydata_nanda),
    "matched":        mydata_nanda["in_nanda"].sum(),
    "match_rate_pct": round(100 * mydata_nanda["in_nanda"].mean(), 1),
})

# The same number, by year. 2023 onward matches because we carried Social
# Services forward; without that step these rows would show 0%.
(mydata_nanda.groupby("year")["in_nanda"]
   .agg(rows="size", matched="sum", match_rate_pct=lambda s: round(100 * s.mean(), 1))
   .reset_index())

# anti_join() keeps the rows of the left table that found no partner: the
# complement of the join.
unmatched <- mydata_nanda %>%
  anti_join(
    socialservices_zcta10 %>% select(zcta10, year),
    by = c("zcta10", "year")
  )

# A row can fail for four different reasons, and each calls for a different
# response: the first two are a data-collection problem, the third is a
# coverage limit, the fourth is usually a boundary mismatch. Three occur in
# this file; the third cannot, because we carried Social Services forward.
# The case stays in the table in case that step is removed.
unmatched %>%
  mutate(reason = case_when(
    is.na(zip)    ~ "no ZIP to start from",
    is.na(zcta10) ~ "ZIP not in the crosswalk",
    year > 2022   ~ "year past NaNDA's coverage",
    TRUE          ~ "ZCTA not in the NaNDA file"
  )) %>%
  tabyl(reason)

reason	n	percent
ZCTA not in the NaNDA file	2	0.0057803
ZIP not in the crosswalk	2	0.0057803
no ZIP to start from	342	0.9884393

* The rows that found no partner: the complement of the join. in_nanda == 0 is
* exactly the "ours only" side of the merge.
use `mydata_nanda', clear
keep if in_nanda == 0

* A row can fail for four different reasons, and each calls for a different
* response: the first two are a data-collection problem, the third is a
* coverage limit, the fourth is usually a boundary mismatch. Three occur in
* this file; the third cannot, because we carried Social Services forward.
* The case stays in the table in case that step is removed.
generate str reason = "ZCTA not in the NaNDA file"
replace  reason = "year past NaNDA's coverage" if year > 2022
replace  reason = "ZIP not in the crosswalk"   if missing(zcta10)
replace  reason = "no ZIP to start from"       if missing(zip)
tabulate reason

# The rows that found no partner: the complement of the join.
# indicator="left_only" is the anti-join.
unmatched = (
    mydata_nanda.merge(
        socialservices_zcta10[["zcta10", "year"]].drop_duplicates(),
        how="left", on=["zcta10", "year"], indicator=True,
    )
    .query('_merge == "left_only"')
    .drop(columns="_merge")
)

# A row can fail for four different reasons, and each calls for a different
# response: the first two are a data-collection problem, the third is a
# coverage limit, the fourth is usually a boundary mismatch. Three occur in
# this file; the third cannot, because we carried Social Services forward.
# The case stays in the table in case that step is removed.
unmatched.assign(reason=np.select(
    [unmatched["zip"].isna(), unmatched["zcta10"].isna(), unmatched["year"] > 2022],
    ["no ZIP to start from", "ZIP not in the crosswalk", "year past NaNDA's coverage"],
    default="ZCTA not in the NaNDA file",
))["reason"].value_counts()

Three reasons appear, and they call for different responses:

No ZIP is the big one and a data-collection problem. Nothing downstream fixes it.
ZIP not in the crosswalk is a typo or a ZIP newer than the crosswalk.
ZCTA not in the NaNDA file usually means the two were drawn on different boundaries.

At two rows apiece the last two vanish inside a match rate. In real data they point to a misalignment upstream, which is why the table lists reasons and not just a count. The fourth reason, year past coverage, is absent only because we carried Social Services forward. Without that step, 407 rows would sit in that line. Either way it is a limit of the NaNDA file, and a decision to document. If the reasons table shows anything beyond what Section 3 already explained, run the who-did-you-lose comparison again on in_nanda.

A toy model

A single model, to show the join produced something usable: does owning a pet go with the neighborhood measures we just attached? The service count is cut into bins so the coefficients read as comparisons against ZCTAs with none. The model drops every unmatched row without a warning, so check the number of rows it used.

mydata_nanda <- mydata_nanda %>%
  mutate(count_totindivfamilyservices_6cat = cut(
    count_totindivfamilyservices,
    breaks = c(-1, 0, 1, 2, 3, 5, 10, Inf),
    labels = c("0", "1", "2", "3", "4-5", "6-10", "11+"),
    right = TRUE
  ))

own_pet_model <- glm(
  own_pet ~ hot_meal_days + count_totindivfamilyservices_6cat + ANY_OPEN_PARK + AFFLUENCE,
  data = mydata_nanda, family = binomial
)
summary(own_pet_model)


Call:
glm(formula = own_pet ~ hot_meal_days + count_totindivfamilyservices_6cat + 
    ANY_OPEN_PARK + AFFLUENCE, family = binomial, data = mydata_nanda)

Coefficients:
                                       Estimate Std. Error z value Pr(>|z|)   
(Intercept)                            0.888723   0.294476   3.018  0.00254 **
hot_meal_days                         -0.076618   0.037718  -2.031  0.04222 * 
count_totindivfamilyservices_6cat1     0.406203   0.256884   1.581  0.11382   
count_totindivfamilyservices_6cat2     0.129286   0.276773   0.467  0.64041   
count_totindivfamilyservices_6cat3     0.736082   0.322711   2.281  0.02255 * 
count_totindivfamilyservices_6cat4-5  -0.008369   0.248372  -0.034  0.97312   
count_totindivfamilyservices_6cat6-10  0.026074   0.228034   0.114  0.90897   
count_totindivfamilyservices_6cat11+   0.084177   0.201700   0.417  0.67643   
ANY_OPEN_PARK                          0.290443   0.195731   1.484  0.13784   
AFFLUENCE                             -0.735152   0.347754  -2.114  0.03452 * 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 2170.5  on 1656  degrees of freedom
Residual deviance: 2150.7  on 1647  degrees of freedom
  (347 observations deleted due to missingness)
AIC: 2170.7

Number of Fisher Scoring iterations: 4

use `mydata_nanda', clear
generate byte count_totindivfamilyservices_6cat = .
replace count_totindivfamilyservices_6cat = 1 if count_totindivfamilyservices == 0
replace count_totindivfamilyservices_6cat = 2 if count_totindivfamilyservices == 1
replace count_totindivfamilyservices_6cat = 3 if count_totindivfamilyservices == 2
replace count_totindivfamilyservices_6cat = 4 if count_totindivfamilyservices == 3
replace count_totindivfamilyservices_6cat = 5 if inrange(count_totindivfamilyservices, 4, 5)
replace count_totindivfamilyservices_6cat = 6 if inrange(count_totindivfamilyservices, 6, 10)
replace count_totindivfamilyservices_6cat = 7 if count_totindivfamilyservices >= 11 & !missing(count_totindivfamilyservices)
label define svc6 1 "0" 2 "1" 3 "2" 4 "3" 5 "4-5" 6 "6-10" 7 "11+"
label values count_totindivfamilyservices_6cat svc6

* ANY_OPEN_PARK is a labeled factor in the .dta; i. treats it as categorical.
logit own_pet hot_meal_days i.count_totindivfamilyservices_6cat i.ANY_OPEN_PARK AFFLUENCE

mydata_nanda["count_totindivfamilyservices_6cat"] = pd.cut(
    mydata_nanda["count_totindivfamilyservices"],
    bins=[-1, 0, 1, 2, 3, 5, 10, np.inf],
    labels=["0", "1", "2", "3", "4-5", "6-10", "11+"],
)

# pandas has no models. statsmodels is the one extra package this step needs
# (pip install statsmodels); everything above runs without it.
try:
    import statsmodels.formula.api as smf
except ImportError:
    smf = None
    print("statsmodels is not installed; skipping the toy model.")

if smf is not None:
    own_pet_model = smf.logit(
        "own_pet ~ hot_meal_days + C(count_totindivfamilyservices_6cat) + C(ANY_OPEN_PARK) + AFFLUENCE",
        data=mydata_nanda.dropna(subset=["count_totindivfamilyservices_6cat", "ANY_OPEN_PARK", "AFFLUENCE"]),
    ).fit()
    print(own_pet_model.summary())

Four questions to ask before every link

What location information do my data contain? An address, a ZIP code, a Census code?
Which census year’s boundaries do my codes follow, and which does the NaNDA file follow?
Which crosswalk did I use, and which boundaries does it produce?
How many rows matched, and do the unmatched rows look different from the rest?

5. Now what?

Grace Noppert closes with a few minutes on why the link is worth making. Your data plus NaNDA data, joined on one geographic identifier, is a file that can answer questions neither could alone. The three measures we attached today each open one:

Social connection. Do older adults near more social services report less loneliness?
Physical activity. Is park access tied to staying active with age?
Food security. Does neighborhood disadvantage shape who goes without a hot meal?

Two cautions:

A neighborhood measure describes a place. Be clear when it is context and when it is a proxy for exposure.
People in the same ZCTA are not independent observations. The analysis has to account for that.

If one of these questions is close to yours, NaNDA has the measure and the catalog is where to start.

The materials ZIP holds:

the synthetic data and its codebook
the crosswalk
this notebook as a runnable notebook.qmd, and the same steps as standalone scripts in r/, stata/ and python/, under the same block headers. Each folder also has a tract version of the merge, for data that already carries a 2010 tract ID
an offline copy of this page

The NaNDA files are the one thing not in the ZIP. They come from ICPSR, which asks for a free account, so you download your own. It takes about five minutes, and data/nanda/README.md in the ZIP walks you through it. You only need them to re-run the code. If you use these materials in teaching or in a paper, cite them by the DOI below. It always resolves to the latest version.

Take this with you

This notebook	Linking NaNDA With Your Data: notebook
Everything, zipped	Workshop materials (ZIP)
Source	GitHub repository
Cite it	DOI 10.5281/zenodo.23071407

Where to go next

The NaNDA catalog on openICPSR. Every dataset, every geography, free to download. When a NaNDA measure ends up in something you publish, cite the dataset itself by its DOI, not just these materials. Each dataset has a ready-made citation on its ICPSR page.
The NaNDA website. Documentation, news, and how to reach the team.
Our published merge code: Chenoweth & Khan, code for merging NaNDA ZCTA-level datasets with the UDS Mapper crosswalk (2024).
The crosswalk: UDS Mapper ZIP Code to ZCTA Crosswalk, 2019 release, archived copy.
Questions after today: nanda-admin@umich.edu

Session information

sessionInfo()

R version 4.4.3 (2025-02-28 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 11 x64 (build 26200)

Matrix products: default


locale:
[1] LC_COLLATE=English_United States.utf8 
[2] LC_CTYPE=English_United States.utf8   
[3] LC_MONETARY=English_United States.utf8
[4] LC_NUMERIC=C                          
[5] LC_TIME=English_United States.utf8    

time zone: America/Denver
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] janitor_2.2.1   here_1.0.1      readxl_1.4.5    lubridate_1.9.4
 [5] forcats_1.0.0   stringr_1.5.1   dplyr_1.1.4     purrr_1.1.0    
 [9] readr_2.1.5     tidyr_1.3.1     tibble_3.3.0    ggplot2_3.5.2  
[13] tidyverse_2.0.0

loaded via a namespace (and not attached):
 [1] bit_4.6.0          gtable_0.3.6       jsonlite_2.0.0     crayon_1.5.3      
 [5] compiler_4.4.3     tidyselect_1.2.1   parallel_4.4.3     snakecase_0.11.1  
 [9] scales_1.4.0       yaml_2.3.10        fastmap_1.2.0      R6_2.6.1          
[13] generics_0.1.4     knitr_1.50         htmlwidgets_1.6.4  rprojroot_2.1.1   
[17] pillar_1.11.0      RColorBrewer_1.1-3 tzdb_0.5.0         rlang_1.1.6       
[21] stringi_1.8.7      xfun_0.53          bit64_4.6.0-1      timechange_0.3.0  
[25] cli_3.6.5          withr_3.0.2        magrittr_2.0.3     digest_0.6.37     
[29] grid_4.4.3         vroom_1.6.5        hms_1.1.3          lifecycle_1.0.4   
[33] vctrs_0.6.5        evaluate_1.0.4     glue_1.8.0         cellranger_1.1.0  
[37] farver_2.1.2       rmarkdown_2.29     tools_4.4.3        pkgconfig_2.0.3   
[41] htmltools_0.5.8.1

about

import sys
print(sys.version)
print("pandas", pd.__version__, "| numpy", np.__version__)