I thought this was an interesting question, so I wrote a little R code to pull PUMS via tidycensus and create these estimates. If you are not familiar, this is a good walk-through of how to use tidycensus for PUMS data. If you also want standard errors with these estimates, you can download the replicate weights as well -- there are examples of how to do this in the linked tidycensus article.
The code below uses 2019 PUMS data, but it looks like 2021 1-year PUMS is set to be released later this week. This code should work exactly the same if you change the year to 2021. Also, I choose to use the ESR (employment status variable) to find people in the household who are employed, but depending on how you want to calculate this, you could use WAGP (person wages) or another measure to identify how many household members have incomes.
library(tidycensus)
library(tidyverse)
# get pums data from api
pums variables = c("AGEP", "ESR"), # person age and employment status
state = "VT",
year = 2019,
survey = "acs1"
)
# aggregate to hh-level
pums_by_hh
filter(!str_detect(SERIALNO, "GQ")) |> # drop group quarters
group_by(SERIALNO, WGTP) |>
filter(any(AGEP <= 6)) |> # keep only hhs with a 6 year old or younger
summarize(
n_tot = n(), # hh size
n_6_under = sum(AGEP <= 6), # number under 6 in hh
n_emp = sum(ESR %in% c("1", "2", "4")) # number of employed people in hh
) |>
ungroup()
# count households by number of people employed
pums_by_hh |>
count(n_emp, wt = WGTP)
# A tibble: 6 × 2
n_emp n
1 0 1033
2 1 9587
3 2 17137
4 3 1220
5 4 66
6 6 186