The information from Spacedman is correct, but Matt is asking how to execute a data join between the douglas_zips feature and a separate table containing a list of lead levels.
The 'join' is an essential function in GIS and relational databases broadly, so it's a critical skill.
We want to take our douglas_zips feature and match all the records from the lead_levels table by the zipcode, adding the columns from lead_levels to douglas_zips. We can use merge() to accomplish that.
https://www.rdocumentation.org/packages/base/versions/3.6.2/topics/merge
NOTE: specify either 'cb=TRUE, year=2020' or 'cb=FALSE' (the default) for ZCTAs from tigris.
library(tigris)
library(sf)
library(plyr)
douglas_zips <- zctas(cb=TRUE, year=2020, starts_with=c("80108","80109", "80104", "80116", "80126", "80129", "80130", "80118", "80124", "80131", "80134", "80138", "80125", "80135"))
## made up data.frame to simulate what you described as your dataset
zipcode = c("80108","80109", "80104", "80116", "80126", "80129", "80130", "80118", "80124", "80131", "80134", "80138", "80125", "80135")
town = c("town1","town2","town3","town4","town5","town6","town7","town8","town9","town10","town11","town12","town13","town14")
lead_level = c(1,2,3,4,5,6,7,8,9,10,11,12,13,14)
dataset = data.frame(zipcode, town, lead_level)
douglas_lead_levels <- merge(douglas_zips, dataset, by.x="ZCTA5CE20",by.y="zipcode")
The result will add two columns to douglas_zips: "town" and "lead_level".
From there, we can plot the lead_level variable as Spacedman explained:
plot(douglas_lead_levels["lead_level"])