if (rebuild_processed) {
# Parse only the variables used for filtering, cohort assignment, weights,
# and fertility counts.
ipums_data <- read_ipums_micro(
ddi,
vars = analysis_vars,
data_file = dat_gz_file,
verbose = FALSE
) %>%
mutate(across(all_of(analysis_vars), as.numeric))
# --- Main series: age-first cohort assignment (membership defined by AGE within
# a census). The 1950 cohort is a proxy from women aged 60-69 in 1960; direct
# 1950 (women 50-59 in 1950) is held out for validation. Birth-year spans are
# derived afterwards for labeling only.
cohort_records <- ipums_data %>%
filter(SEX == 2) %>%
mutate(
source_year = YEAR,
cohort_source = case_when(
YEAR == 1950 & between(AGE, 50, 59) ~ "direct_1950_validation",
YEAR == 1960 & between(AGE, 60, 69) ~ "proxy_1950",
YEAR == 1960 & between(AGE, 50, 59) ~ "main_1960",
YEAR == 1970 & between(AGE, 50, 59) ~ "main_1970",
YEAR == 1980 & between(AGE, 50, 59) ~ "main_1980",
YEAR == 1990 & between(AGE, 50, 59) ~ "main_1990",
TRUE ~ NA_character_
),
Year = case_when(
cohort_source %in% c("proxy_1950", "direct_1950_validation") ~ 1950,
cohort_source == "main_1960" ~ 1960,
cohort_source == "main_1970" ~ 1970,
cohort_source == "main_1980" ~ 1980,
cohort_source == "main_1990" ~ 1990,
TRUE ~ NA_real_
),
series = case_when(
cohort_source == "direct_1950_validation" ~ "validation",
!is.na(cohort_source) ~ "main",
TRUE ~ NA_character_
),
valid_chborn = CHBORN > 0 & CHBORN < 97,
n_child = if_else(valid_chborn, pmin(CHBORN - 1, 12), NA_real_)
) %>%
filter(!is.na(cohort_source))
cohort_composition <- cohort_records %>%
group_by(series, cohort_source, source_year, Year) %>%
summarise(
n_records = n(),
age_min = min(AGE, na.rm = TRUE),
age_max = max(AGE, na.rm = TRUE),
birthyr_min = min(BIRTHYR[BIRTHYR > 1800], na.rm = TRUE),
birthyr_max = max(BIRTHYR[BIRTHYR > 1800], na.rm = TRUE),
.groups = "drop"
) %>%
mutate(cohort = paste0(birthyr_min, "-", birthyr_max)) %>%
arrange(Year, source_year, cohort_source)
cohort_records <- cohort_records %>%
left_join(select(cohort_composition, cohort_source, cohort), by = "cohort_source")
age_cohort_validation <- cohort_composition
completed_fertility <- cohort_records
valid_completed_fertility <- completed_fertility %>% filter(valid_chborn)
main_fertility <- valid_completed_fertility %>% filter(series == "main")
validation_1950 <- valid_completed_fertility %>%
filter(cohort_source %in% c("proxy_1950", "direct_1950_validation"))
# How much data the CHBORN validity rule excludes, by cohort source.
chborn_universe <- completed_fertility %>%
count(cohort_source, source_year, Year, cohort, valid_chborn, name = "n_records") %>%
pivot_wider(
names_from = valid_chborn,
values_from = n_records,
names_prefix = "valid_chborn_",
values_fill = 0
) %>%
mutate(
total_records = valid_chborn_FALSE + valid_chborn_TRUE,
valid_share = valid_chborn_TRUE / total_records
) %>%
arrange(Year, source_year, cohort_source)
cohort_means <- valid_completed_fertility %>%
group_by(cohort_source, series, source_year, Year, cohort) %>%
summarise(
n_women = n(),
n_women_weighted = sum(PERWT, na.rm = TRUE),
mean_children = mean(n_child),
mean_children_weighted = weighted.mean(n_child, PERWT, na.rm = TRUE),
share_12plus = mean(n_child == 12),
.groups = "drop"
) %>%
arrange(series, Year, source_year, cohort_source)
fertility_counts <- main_fertility %>%
count(Year, cohort, n_child, name = "n_women") %>%
arrange(Year, n_child)
fertility_counts_weighted <- main_fertility %>%
group_by(Year, cohort, n_child) %>%
summarise(n_women_weighted = sum(PERWT, na.rm = TRUE), .groups = "drop") %>%
arrange(Year, n_child)
fertility_pmf <- fertility_counts %>%
group_by(Year, cohort) %>%
mutate(probability = n_women / sum(n_women)) %>%
ungroup() %>%
transmute(Year, cohort, children = n_child, probability)
fertility_pmf_weighted <- fertility_counts_weighted %>%
group_by(Year, cohort) %>%
mutate(probability_weighted = n_women_weighted / sum(n_women_weighted)) %>%
ungroup() %>%
transmute(Year, cohort, children = n_child, probability_weighted)
cohort_sample_sizes <- main_fertility %>%
group_by(Year, cohort) %>%
summarise(
n_women = n(),
n_women_weighted = sum(PERWT, na.rm = TRUE),
.groups = "drop"
) %>%
arrange(Year)
validation_1950_summary <- validation_1950 %>%
group_by(cohort_source) %>%
summarise(
source_year = first(source_year),
output_year = first(Year),
n_women = n(),
n_women_weighted = sum(PERWT, na.rm = TRUE),
mean_children = mean(n_child),
mean_children_weighted = weighted.mean(n_child, PERWT, na.rm = TRUE),
share_12plus = mean(n_child == 12),
.groups = "drop"
) %>%
arrange(cohort_source)
validation_1950_pmf <- validation_1950 %>%
count(cohort_source, n_child, name = "n_women") %>%
left_join(
validation_1950 %>%
group_by(cohort_source, n_child) %>%
summarise(n_women_weighted = sum(PERWT, na.rm = TRUE), .groups = "drop"),
by = c("cohort_source", "n_child")
) %>%
group_by(cohort_source) %>%
mutate(
probability = n_women / sum(n_women),
probability_weighted = n_women_weighted / sum(n_women_weighted)
) %>%
ungroup() %>%
arrange(cohort_source, n_child)
# --- Common-universe analysis (separate frame; the main series is untouched).
# Follow decade birth cohorts across the 1960-1990 censuses and, where the whole
# cohort has reached age 50+, contrast the children-ever-born distribution under
# two universes: all women vs ever-married. This isolates the effect of the
# universe shift (1950/1960 asked only ever-married women; 1970-1990 all women).
# MARST: 1-5 = ever-married, 6 = never-married, 9 = blank/missing (kept in the
# all-women universe but excluded from the ever-married subset).
band_records <- ipums_data %>%
filter(SEX == 2, between(YEAR, 1960, 1990), between(BIRTHYR, 1890, 1939)) %>%
mutate(
band_start = 10L * (BIRTHYR %/% 10L),
birth_band = paste0(band_start, "-", band_start + 9L),
universe_year = YEAR,
# Keep only cells where the whole decade cohort has reached age 50+.
cohort_complete = (universe_year - (band_start + 9L)) >= 50L,
valid_chborn = CHBORN > 0 & CHBORN < 97,
n_child = if_else(valid_chborn, pmin(CHBORN - 1, 12), NA_real_),
marital = case_when(
MARST %in% 1:5 ~ "ever_married",
MARST == 6 ~ "never_married",
TRUE ~ "unknown"
),
gq_household = GQ %in% c(1, 2, 5)
) %>%
filter(cohort_complete, valid_chborn)
band_long <- bind_rows(
band_records %>% mutate(universe = "all_women"),
band_records %>% filter(marital == "ever_married") %>% mutate(universe = "ever_married")
)
band_summary <- band_long %>%
group_by(birth_band, universe_year, universe) %>%
summarise(
age_min = min(AGE),
age_max = max(AGE),
n_women = n(),
n_women_weighted = sum(PERWT),
mean_children = mean(n_child),
mean_children_weighted = weighted.mean(n_child, PERWT),
var_children = var(n_child),
share_0 = mean(n_child == 0),
share_12plus = mean(n_child == 12),
gq_household_share = mean(gq_household),
.groups = "drop"
) %>%
arrange(birth_band, universe_year, universe)
band_pmf <- band_long %>%
count(birth_band, universe_year, universe, n_child, name = "n_women") %>%
left_join(
band_long %>%
group_by(birth_band, universe_year, universe, n_child) %>%
summarise(n_women_weighted = sum(PERWT), .groups = "drop"),
by = c("birth_band", "universe_year", "universe", "n_child")
) %>%
group_by(birth_band, universe_year, universe) %>%
mutate(
probability = n_women / sum(n_women),
probability_weighted = n_women_weighted / sum(n_women_weighted)
) %>%
ungroup() %>%
arrange(birth_band, universe_year, universe, n_child)
# Within-census contrast. The two universes coincide by construction before 1970
# (never-married women were not in the CHBORN universe), so the contrast is
# meaningful only for universe_year >= 1970.
universe_contrast <- band_summary %>%
select(birth_band, universe_year, universe, n_women, mean_children, var_children, share_0) %>%
pivot_wider(
names_from = universe,
values_from = c(n_women, mean_children, var_children, share_0)
) %>%
mutate(
mean_shift = mean_children_all_women - mean_children_ever_married,
share0_shift = share_0_all_women - share_0_ever_married,
pct_never_married = 1 - n_women_ever_married / n_women_all_women
)
tvd_table <- band_pmf %>%
select(birth_band, universe_year, universe, n_child, probability, probability_weighted) %>%
pivot_wider(
names_from = universe,
values_from = c(probability, probability_weighted),
values_fill = 0
) %>%
group_by(birth_band, universe_year) %>%
summarise(
tvd = 0.5 * sum(abs(probability_all_women - probability_ever_married)),
tvd_weighted = 0.5 * sum(abs(probability_weighted_all_women - probability_weighted_ever_married)),
.groups = "drop"
)
universe_contrast <- universe_contrast %>%
left_join(tvd_table, by = c("birth_band", "universe_year")) %>%
arrange(birth_band, universe_year)
# Early-cohort decomposition: the two cohorts seen ever-married at completed
# fertility (1960) and all-women later (1970). Split the gap into the universe
# effect (same census, adding never-married women) and the mortality/decade
# effect (same ever-married universe, older age).
pick_mean <- function(b, y, u) {
v <- band_summary$mean_children[
band_summary$birth_band == b &
band_summary$universe_year == y &
band_summary$universe == u
]
if (length(v) == 1) v else NA_real_
}
early_cohort_followup <- tibble(
birth_band = c("1890-1899", "1900-1909"),
completed_year = c(1960, 1960),
later_year = c(1970, 1970)
) %>%
mutate(
completed_em_mean = mapply(function(b, y) pick_mean(b, y, "ever_married"), birth_band, completed_year),
later_em_mean = mapply(function(b, y) pick_mean(b, y, "ever_married"), birth_band, later_year),
later_all_mean = mapply(function(b, y) pick_mean(b, y, "all_women"), birth_band, later_year),
universe_effect = later_all_mean - later_em_mean,
mortality_effect = later_em_mean - completed_em_mean
)
# Self-check: the 1900-1909 / 1960 ever-married mean should be close to the main
# 1960 series mean, and 1890-1899 / 1960 close to the proxy-1950 mean.
message(
"Self-check (common-universe vs main series): ",
"1900-1909/1960 ever-married mean = ", round(pick_mean("1900-1909", 1960, "ever_married"), 3),
" (main_1960 ~ 2.429); ",
"1890-1899/1960 ever-married mean = ", round(pick_mean("1890-1899", 1960, "ever_married"), 3),
" (proxy_1950 ~ 2.818)."
)
# Write the aggregate CSVs.
write_csv(cohort_sample_sizes, file.path(processed_dir, "ipums_cohort_sample_sizes.csv"))
write_csv(fertility_counts, file.path(processed_dir, "ipums_fertility_counts_by_cohort.csv"))
write_csv(fertility_pmf, file.path(processed_dir, "ipums_fertility_pmf_by_cohort.csv"))
write_csv(fertility_counts_weighted, file.path(processed_dir, "ipums_fertility_counts_by_cohort_weighted.csv"))
write_csv(fertility_pmf_weighted, file.path(processed_dir, "ipums_fertility_pmf_by_cohort_weighted.csv"))
write_csv(validation_1950_summary, file.path(processed_dir, "validation_1950_direct_vs_proxy.csv"))
write_csv(validation_1950_pmf, file.path(processed_dir, "validation_1950_direct_vs_proxy_pmf.csv"))
write_csv(age_cohort_validation, file.path(processed_dir, "ipums_age_cohort_validation.csv"))
write_csv(chborn_universe, file.path(processed_dir, "ipums_chborn_universe_by_cohort.csv"))
write_csv(cohort_means, file.path(processed_dir, "ipums_cohort_fertility_checks.csv"))
write_csv(band_summary, file.path(processed_dir, "ipums_cohort_universe_summary.csv"))
write_csv(band_pmf, file.path(processed_dir, "ipums_cohort_universe_pmf.csv"))
write_csv(universe_contrast, file.path(processed_dir, "ipums_universe_contrast_by_cohort.csv"))
write_csv(early_cohort_followup, file.path(processed_dir, "ipums_early_cohort_followup.csv"))
write_csv(
tibble(
processing_version = processing_version,
cohort_definition = "main series: age-defined cohorts (50-59; proxy-1950 from 60-69 in 1960; direct-1950 validation-only); CHBORN harmonized to 12+. Common-universe analysis: decade birth cohorts followed across 1960-1990 at completed ages, all-women vs ever-married (MARST 1-5).",
rebuilt_at = as.character(Sys.time())
),
file.path(processed_dir, "ipums_preprocess_metadata.csv")
)
fertility_pmf %>%
group_split(Year) %>%
lapply(function(dat) {
year <- unique(dat$Year)
write_csv(
dat %>% transmute(children, probability),
file.path(processed_dir, paste0("fertility_pmf_", year, ".csv"))
)
})
rm(
ipums_data, cohort_records, cohort_composition, completed_fertility,
valid_completed_fertility, main_fertility, validation_1950,
band_records, band_long, tvd_table
)
} else {
# Load the committed aggregates without reparsing restricted person-level data.
rd <- function(f) read_csv(file.path(processed_dir, f), show_col_types = FALSE)
cohort_sample_sizes <- rd("ipums_cohort_sample_sizes.csv")
fertility_counts <- rd("ipums_fertility_counts_by_cohort.csv")
fertility_pmf <- rd("ipums_fertility_pmf_by_cohort.csv")
fertility_counts_weighted <- rd("ipums_fertility_counts_by_cohort_weighted.csv")
fertility_pmf_weighted <- rd("ipums_fertility_pmf_by_cohort_weighted.csv")
validation_1950_summary <- rd("validation_1950_direct_vs_proxy.csv")
validation_1950_pmf <- rd("validation_1950_direct_vs_proxy_pmf.csv")
age_cohort_validation <- rd("ipums_age_cohort_validation.csv")
chborn_universe <- rd("ipums_chborn_universe_by_cohort.csv")
cohort_means <- rd("ipums_cohort_fertility_checks.csv")
band_summary <- rd("ipums_cohort_universe_summary.csv")
band_pmf <- rd("ipums_cohort_universe_pmf.csv")
universe_contrast <- rd("ipums_universe_contrast_by_cohort.csv")
early_cohort_followup <- rd("ipums_early_cohort_followup.csv")
preprocess_metadata <- rd("ipums_preprocess_metadata.csv")
if (!identical(preprocess_metadata$processing_version[[1]], processing_version)) {
stop(
"Processed files were built with processing version ",
preprocess_metadata$processing_version[[1]],
"; expected ", processing_version,
". Re-render with KINCOUNTS_REBUILD_IPUMS=true.",
call. = FALSE
)
}
}
# Confirm the complete aggregate output set is available after either branch.
written_files <- tibble(
file = processed_files,
status = if_else(file.exists(processed_paths), "available", "missing")
)
if (any(written_files$status == "missing")) {
stop(
"Missing processed files: ",
paste(written_files$file[written_files$status == "missing"], collapse = ", "),
call. = FALSE
)
}
