Skip to content

sync with main #19

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
kyleGrealis wants to merge 2 commits into from
Closed

sync with main #19

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
7e25ca4
Harmonize ALQ across cycles; correct PAQ 2017 coverage; add pad680.
Goulett May 9, 2026
6d40725
Fix ALQ and OCQ cross-cycle harmonization in create_base_expanded.R:
Goulett May 11, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
162 changes: 125 additions & 37 deletions inst/scripts/derived/create_base_expanded.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -57,11 +57,26 @@
# TCHOL – Available 2005+ only (earlier cycles used l13/l13_b/l13_c).
# PAQ – The paq605/pad615 variable series used here was introduced in
# 2007–2008. Cycles 1999–2006 used different variable names not
# harmonized here, and the series was removed again in 2017–2018
# and 2021. total_met_min is NA by design for year==1999, 2001,
# 2003, 2005, 2017, and 2021. Valid data: 2007–2015 only.
# ALQ – Not collected 1999–2000; adults 20+ only → high structural
# missingness in the full sample.
# harmonized here, so total_met_min is NA by design for those
# cycles. 2021 uses an entirely different schema (pad790q/u,
# pad810q/u sedentary-time module) — also NA by design.
# 2017–2018 uses the SAME paq605–pad675 schema as 2007–2015
# (gates at 100%, sub-questions at 24–42%) — total_met_min IS
# computed for 2017. Valid data: 2007–2017 (6 cycles).
# pad680 (sedentary minutes/day) is retained separately; it is
# available for 2007–2021 and not used in the MET-min formula.
# ALQ – Collected all cycles 1999–2017; adults 20+ only → ~6–13%
# structural missingness per cycle (age restrictions, refusals).
# No full-cycle design gap exists.
# alq110 ("Had ≥12 drinks in any one year?") was renamed to
# alq111 in 2017–2018 and 2021. alq110 is NA for those cycles;
# alq111 is NA for all earlier cycles. Both are coalesced into
# alq110_harmonized in Section 4 so alcohol_status is derived
# consistently across all cycles.
# LIMITATION: the frequency follow-up alq120q was also
# restructured in 2017 (→ alq121); alq121 is not harmonized
# here. 2017/2021 drinkers (alq111 == "Yes") will therefore be
# coded "Not provided" rather than "Current" or "Former".
# DUQ – Not available for 1999, 2001, or 2003 (NA by design).
# 1999: not collected. 2001 (DUQ_B) and 2003 (DUQ_C): collected
# but NCHS-restricted (sensitive data); absent from nhanesdata.
Expand Down Expand Up @@ -93,16 +108,24 @@
#
# HARMONIZATION NOTES
# -------------------
# Several DEMO variables changed names or response-option wording across cycles.
# Several variables changed names or response-option wording across cycles.
# Key issues handled here:
# income – indhhinc (1999–2006) vs indhhin2 (2007+): coalesced into one
# categorical variable; indfmpir (PIR) used as continuous income.
# income – indhhinc (1999–2006) vs indhhin2 (2007+): coalesced into one
# categorical variable; indfmpir (PIR) used as continuous income.
# nativity – dmdborn/dmdborn2/dmdborn4 changed labels across cycles; we use
# the cleaned nhanesdata label strings (case_when on text).
# race – ridreth3 (Asian category) only exists 2011+; ridreth1 used as
# fallback.
# race – ridreth3 (Asian category) only exists 2011+; ridreth1 used as
# fallback.
# BPX/BPXO – auscultatory (bpxsy*/bpxdi*) coalesced with oscillometric
# (bpxosy*/bpxodi*) to handle the 2021 measurement transition.
# ALQ – alq110 ("Had ≥12 drinks in any one year?") renamed to alq111 in
# 2017–2018 and 2021. Harmonized via:
# alq110_harmonized = coalesce(alq110, alq111)
# This intermediate column drives alcohol_status derivation; it is
# not retained in base_small. alq111 labels confirmed: "Yes"/"No"
# in 2017; "Yes"/"No"/"Don't know" in 2021. The frequency follow-
# up (alq120q → alq121 in 2017+) is NOT harmonized; 2017/2021
# drinkers ("Yes") fall into "Not provided" in alcohol_status.
################################################################################


Expand Down Expand Up @@ -398,19 +421,58 @@ base <- base |>

# -------------------------------------------------------------------------
# ALCOHOL USE
# -------------------------------------------------------------------------
#
# CROSS-CYCLE HARMONIZATION — three distinct schemas:
#
# 1999: ALQ collected; same alq110 gate as 2001–2003.
#
# 2001–2003: alq110 = "In any one year, had ≥12 drinks?" (Yes/No gate).
# alq120q = frequency in past 12 months (numeric, 0–365; 0 →
# former drinker; 1–700 → current; 999 → don't know).
#
# 2005 (ALQ_D): TWO-GATE STRUCTURE — different from all other cycles.
# alq101 = "Ever had ≥12 drinks in any one year?" (primary gate).
# Yes → skip alq110; use alq120q for Current/Former (same scale as 2001–2003).
# No → alq110 = "Had any drink in past 12 months?" (secondary gate).
# alq110 = "Yes" → current light drinker → "Current".
# alq110 = "No" → abstainer → "Never".
# 2005 branches are evaluated FIRST in case_when to prevent them from
# falling through to the generic alq110_harmonized logic below.
#
# 2007–2015: Same alq110 gate as 2001–2003; alq120q same scale.
#
# 2017, 2021: alq110 renamed to alq111 (labels identical: Yes/No/Don't know).
# alq110_harmonized = coalesce(alq110, alq111).
# alq120q restructured in 2017 (→ alq121, not harmonized here);
# 2017/2021 "Yes" respondents fall into "Not provided".
# -------------------------------------------------------------------------
alq110_harmonized = dplyr::coalesce(
as.character(alq110),
as.character(alq111)
),
alcohol_status = factor(
case_when(
alq110 %in% c("Don't know", "Refused") |
alq120q >= 777 ~ "Not provided",
alq110 == "No" ~ "Never",
alq110 == "Yes" & alq120q == 0 ~ "Former",
alq110 == "Yes" & between(alq120q, 1, 700) ~ "Current",
alq110 == "Yes" & (alq120q == 999 | is.na(alq120q)) ~ "Not provided",
is.na(alq110) & !is.na(alq120u) ~ "Current",
is.na(alq110) & alq120q == 0 ~
# ---- 2005: alq101 primary gate; alq110 = past-12-months drinker ----
year == 2005 & alq101 %in% c("Don't know", "Refused") ~ "Not provided",
year == 2005 & alq101 == "Yes" & alq120q == 0 ~ "Former",
year == 2005 & alq101 == "Yes" & dplyr::between(alq120q, 1, 700) ~ "Current",
year == 2005 & alq101 == "Yes" & (alq120q >= 777 | is.na(alq120q)) ~ "Not provided",
year == 2005 & alq101 == "No" & alq110 == "Yes" ~ "Current",
year == 2005 & alq101 == "No" & alq110 == "No" ~ "Never",
year == 2005 & alq101 == "No" & is.na(alq110) ~ "Does not currently drink, hx unknown",
year == 2005 & is.na(alq101) ~ NA_character_,
# ---- All other cycles: alq110_harmonized gate ----
alq110_harmonized %in% c("Don't know", "Refused") |
alq120q >= 777 ~ "Not provided",
alq110_harmonized == "No" ~ "Never",
alq110_harmonized == "Yes" & alq120q == 0 ~ "Former",
alq110_harmonized == "Yes" & between(alq120q, 1, 700) ~ "Current",
alq110_harmonized == "Yes" & (alq120q == 999 |
is.na(alq120q)) ~ "Not provided",
is.na(alq110_harmonized) & !is.na(alq120u) ~ "Current",
is.na(alq110_harmonized) & alq120q == 0 ~
"Does not currently drink, hx unknown",
.default = alq110
.default = alq110_harmonized
),
levels = c(
"Never", "Does not currently drink, hx unknown",
Expand Down Expand Up @@ -608,11 +670,15 @@ base <- base |>
total_cholesterol = lbxtc,

# -------------------------------------------------------------------------
# OCQ SENTINEL CODE RECODE
# ocq180 uses 77777 (refused) and 99999 (don't know) as skip codes.
# These must be set to NA before analysis; they are not real hour values.
# OCQ SENTINEL CODE RECODE + 2001 VARIABLE NAME HARMONIZATION
# In 2001–2002 (OCQ_B), hours/week was stored as ocd180 (OCD prefix)
# rather than ocq180 (OCQ prefix used in all other cycles). Coalesce
# first so 2001 workers are not silently dropped.
# Sentinel codes 77777 (refused) and 99999 (don't know) present in both
# naming variants; set to NA before analysis.
# -------------------------------------------------------------------------
ocq180 = if_else(ocq180 %in% c(77777L, 99999L), NA_real_, as.numeric(ocq180))
ocq180 = dplyr::coalesce(as.numeric(ocq180), as.numeric(ocd180)),
ocq180 = if_else(ocq180 %in% c(77777, 99999), NA_real_, ocq180)

)

Expand Down Expand Up @@ -660,9 +726,11 @@ base_small <- base |>
smoking_status,
alcohol_status, alcohol_current,
total_met_min, pa_level, phys_active,
# NOTE: total_met_min is NA by design for year==1999, 2001, 2003, 2005,
# 2017, and 2021. Valid MET-min data exists for 2007–2015 only
# (5 of 11 cycles).
# NOTE: total_met_min is NA by design for year==1999, 2001, 2003, 2005
# (pre-2007 PAQ schema) and 2021 (different schema). Valid data
# exists for 2007–2017 (6 of 11 cycles); 2017 uses the same
# paq605–pad675 variables as 2007–2015.
pad680, # sedentary time (minutes/day sitting); available 2007–2021

# ----- Body composition ---------------------------------------------------
bmi, # continuous (kg/m²)
Expand Down Expand Up @@ -807,17 +875,37 @@ message("\ndiabetes distribution by year:")
print(base_small |> count(diabetes, year) |>
tidyr::pivot_wider(names_from = year, values_from = n))

# 6e. PA missingness by year
# Expect 100% missing for year==1999, 2001, 2003, 2005 (pre-2007 PAQ names),
# year==2017 and year==2021 (PAQ restructured). Valid data: 2007–2015 only.
# 6e. Alcohol coverage by year
# Expect ~87–95% non-missing across all cycles (ALQ collected 1999–2017).
# Expect "Never" to be populated for 2017 and 2021 (alq111 coalesced).
# Expect 2017/2021 drinkers to appear in "Not provided" (alq120q not harmonized).
# Expect 2005 to have a real distribution (alq101 two-gate fix applied).
message("\nAlcohol status non-missing rate by year:")
print(
base_small |>
group_by(year) |>
summarise(pct_alc_present = round(mean(!is.na(alcohol_status)) * 100, 1), n = n())
)

message("\nalcohol_status distribution for 2017 and 2021:")
print(
base_small |>
filter(year %in% c(2017, 2021)) |>
count(year, alcohol_status)
)

# 6f. PA missingness by year
# Expect 100% missing for year==1999, 2001, 2003, 2005 (pre-2007 PAQ names)
# and 2021 (different schema). 2017 uses the same paq605–pad675 series as
# 2007–2015 and should show non-zero total_met_min. Valid data: 2007–2017.
message("\nPA (total_met_min) missing rate by year:")
print(
base_small |>
group_by(year) |>
summarise(pct_pa_missing = round(mean(is.na(total_met_min)) * 100, 1), n = n())
)

# 6f. Mortality coverage by year
# 6g. Mortality coverage by year
# Expect non-NA mortstat for 1999–2017 eligible adults; 0% for year==2021.
message("\nMortality coverage by year (% with non-NA mortstat):")
print(
Expand All @@ -826,7 +914,7 @@ print(
summarise(pct_linked = round(mean(!is.na(mortstat)) * 100, 1), n = n())
)

# 6g. DUQ coverage by year
# 6h. DUQ coverage by year
# Expect 0% present for year==1999; age-restriction pattern visible in others.
message("\nDUQ (duq200) non-missing rate by year:")
print(
Expand All @@ -835,7 +923,7 @@ print(
summarise(pct_duq_present = round(mean(!is.na(duq200)) * 100, 1), n = n())
)

# 6h. OCQ work schedule (ocq670) coverage by year
# 6i. OCQ work schedule (ocq670) coverage by year
# Expect 0% present for years 1999–2011 (variable collected 2013+ only).
message("\nOCQ work schedule (ocq670) non-missing rate by year:")
print(
Expand All @@ -844,7 +932,7 @@ print(
summarise(pct_ocq670_present = round(mean(!is.na(ocq670)) * 100, 1), n = n())
)

# 6i. Pregnancy variable (ridexprg) coverage by year
# 6j. Pregnancy variable (ridexprg) coverage by year
# Expect NA for all males and age-ineligible participants; non-NA for
# examined females. Overall missingness ~81% in the full sample is expected.
message("\nPregnancy status (ridexprg) non-missing rate by year:")
Expand All @@ -854,19 +942,19 @@ print(
summarise(pct_ridexprg_present = round(mean(!is.na(ridexprg)) * 100, 1), n = n())
)

# 6j. Survey weight sanity checks
# 6k. Survey weight sanity checks
stopifnot(
"Negative wtmec2yr found" = all(base_small$wtmec2yr >= 0, na.rm = TRUE),
"Negative wtint2yr found" = all(base_small$wtint2yr >= 0, na.rm = TRUE)
)

# 6k. PIR range check
# 6l. PIR range check
if (any(base_small$pir > 10, na.rm = TRUE)) {
warning(sum(base_small$pir > 10, na.rm = TRUE),
" rows have PIR > 10 — check for coding errors.")
}

# 6l. OCQ hours check (>99 hrs/week; sentinel codes 77777/99999 already NA).
# 6m. OCQ hours check (>99 hrs/week; sentinel codes 77777/99999 already NA).
# Remaining values >99 are expected to be real extreme-hour workers.
# Note: some cycles used 100 as a ceiling code ("100 or more") — interpret
# the value 100 with caution in continuous models.
Expand Down
Loading