NHANES III (1988-1994): Fixed-Width Files, Cotinine, and Mortality
Source:vignettes/nhanes-iii.Rmd
nhanes-iii.RmdOverview
NHANES III (1988-1994) predates the continuous survey. Its public-use
data are not per-cycle XPT files but a few large fixed-width
.dat files, each with a SAS program (.sas)
giving the column positions. nhanesR handles the NHANES III
mortality file directly
(nhanes_mortality_parse(cycles = "1988-1994")); this
vignette shows how to read the survey files themselves and link
them.
The example builds an adult analysis file with measured height, weight and waist, cigarette smoking status checked against serum cotinine, and mortality follow-up through 2019. NHANES III is worth the extra work for survival analysis: participants have a median of about 27 years of follow-up.
None of the code chunks run automatically. The exam file alone is about 195 MB.
A complete, runnable version of this workflow (download, parse,
recode, link mortality, save an .rds) ships with the
package:
system.file("scripts", "assemble_nhanes3_data.R", package = "nhanesR").
Copy it to a working directory and run it there; it caches the raw files
in ./nhanes3_raw/ and writes
nhanes3_pooled.rds.
1. Download the files
The main files are in nhanes3/1a/. Serum cotinine is in
the second laboratory file, nhanes3/2a/lab2.dat.
The main lab.dat has a COP column, but it is
blank for every record.
raw <- file.path(tempdir(), "nhanes3")
dir.create(raw, showWarnings = FALSE)
options(timeout = 900) # default 60 s is too short
base <- "https://wwwn.cdc.gov/nchs/data/nhanes3/"
files <- c("1a/exam.dat", "1a/exam.sas", "1a/exam-acc.pdf",
"1a/adult.dat", "1a/adult.sas", "1a/adult-acc.pdf",
"2a/lab2.dat", "2a/lab2.sas", "2a/lab2-acc.pdf")
for (f in files) {
dest <- file.path(raw, basename(f))
if (!file.exists(dest)) download.file(paste0(base, f), dest, mode = "wb")
}Each data file has a codebook alongside it, named
<file>-acc.pdf (for example
1a/exam-acc.pdf). The codebooks list value labels and
missing-value codes; section 5 shows how to pull entries out of
them.
2. Read column positions from the SAS program
Each .sas file has an INPUT block of lines
like SEQN 1-5 or HSSEX 15. The helper below
parses that block and checks it against the record length
(LRECL, which counts 2 bytes for the line ending).
SAScii::parse.SAScii() produces the same layout.
sas_layout <- function(sas_file) {
s <- readLines(sas_file, warn = FALSE)
lrecl <- as.integer(sub(".*LRECL\\s*=\\s*([0-9]+).*", "\\1",
grep("LRECL\\s*=", s, value = TRUE)[1]))
s <- s[(grep("^\\s*INPUT\\s*$", s)[1] + 1):length(s)]
s <- s[seq_len(grep("^\\s*;\\s*$", s)[1] - 1)]
m <- regmatches(s, regexec(
"^\\s*([A-Z][A-Z0-9_]*)\\s+\\$?\\s*([0-9]+)(-([0-9]+))?\\s*$", s))
m <- m[lengths(m) == 5]
L <- data.frame(var = vapply(m, `[`, "", 2),
start = as.integer(vapply(m, `[`, "", 3)))
end <- vapply(m, `[`, "", 5)
L$end <- ifelse(nzchar(end), as.integer(end), L$start)
stopifnot(max(L$end) == lrecl - 2)
L
}3. Read only the columns you need
The exam file has more than 2,000 columns. Reading a handful by position is fast; reading all of them is not.
read_nh3 <- function(stem, vars) {
L <- sas_layout(file.path(raw, paste0(stem, ".sas")))
stopifnot(all(vars %in% L$var))
L <- L[match(vars, L$var), ]
x <- readr::read_fwf(file.path(raw, paste0(stem, ".dat")),
readr::fwf_positions(L$start, L$end, L$var),
col_types = readr::cols(.default = "d"))
attr(x, "width") <- setNames(L$end - L$start + 1, L$var)
x
}
exam <- read_nh3("exam", c("SEQN", "HSSEX", "HSAGEIR", "DMARETHN",
"BMPHT", "BMPWT", "BMPBMI", "BMPWAIST",
"WTPFEX6", "SDPPSU6", "SDPSTRA6"))
adult <- read_nh3("adult", c("SEQN", "HAR1", "HAR3", "HAR16", "HAR24", "HAR27"))
lab2 <- read_nh3("lab2", c("SEQN", "COP"))Use the 6-year weights (WTPFEX6, SDPPSU6,
SDPSTRA6) for the whole survey. The *1 and
*2 versions cover phase 1 (1988-1991) or phase 2
(1991-1994) alone.
The same steps with the package functions
The helper code above is what nhanesR provides as
functions, with variable labels and a search added. You can stop after
the layout, which needs only the small .sas program, and
decide which columns to read before you download the large
.dat file:
nhanes3_files() # the files and their CDC addresses
nhanes3_layout("adult", pattern = "tall|weigh") # variables whose name or label matches
nhanes3_layout("exam", pattern = "^BMP") # the body-measure variables, with labels
x <- nhanes3_read("exam", c("HSSEX", "HSAGEIR", "BMPHT", "BMPWT", "WTPFEX6"),
recode = c("BMPHT", "BMPWT")) # SEQN is always included
attr(x, "var_labels") # the SAS labels travel with the dataIn an interactive session, leave vars out and
nhanes3_read() offers the (optionally
pattern-narrowed) list to pick from. recode
applies only to the variables you name, because IDs, weights and age
have no all-8s/all-9s codes and many items have special codes of their
own (section 5).
4. Missing-value codes
NHANES III stores “blank but applicable” as all 8s and “don’t know”
as all 9s, sized to the field: 8/9 for a
one-character item, 88888 for a five-character measurement.
Left in place they look like real values: an 88888 cotinine code makes
the mean serum cotinine about 10,000 ng/mL.
recode_codes <- function(x, vars) {
w <- attr(x, "width")
for (v in vars) {
codes <- as.numeric(c(strrep("8", w[[v]]), strrep("9", w[[v]])))
x[[v]][x[[v]] %in% codes] <- NA
}
x
}
exam <- recode_codes(exam, c("BMPHT", "BMPWT", "BMPBMI", "BMPWAIST"))
adult <- recode_codes(adult, c("HAR1", "HAR3", "HAR16", "HAR24", "HAR27"))
lab2 <- recode_codes(lab2, "COP")Check each variable’s codebook entry before applying this rule to it
(next section). Age (HSAGEIR) is top-coded at 90 and has no
such codes.
Missing serum cotinine is mostly missing labs, not a missing assay: children under 4 were not eligible, and most people coded 88888 have no serum results at all.
5. Look up variables in the codebook
The all-8s/all-9s rule is not the whole story. Many items have their
own special codes, and their width varies: cigarettes per day
(HAR4S) uses 666 for “varies” and 777 for “less than 1 per
day”; cigars per day (HAR25) uses 77; cigarettes per day at
the heaviest (HAR7S) uses 6666; age started smoking
(HAR2) uses 000 for “never”. Left in place, a 777 reads as
777 cigarettes a day.
The codebook PDFs are plain fixed-width text, one entry per variable:
1336 Is your home drinking water bottled or
HFE4 from the tap (faucet)?
2443 1 Bottled (HFE7)
88 8 Blank but applicableThe helpers below (they need the pdftools package)
return each requested variable’s file position, question text,
interviewer notes, and a table of codes, labels and counts. Counts are
for the whole file, so they also show how many records each special code
affects.
if (!requireNamespace("pdftools", quietly = TRUE)) install.packages("pdftools")
library(pdftools)
nh3_codebook_lines <- function(pdf) {
unlist(lapply(pdftools::pdf_text(pdf), function(p) {
L <- strsplit(p, "\n", fixed = TRUE)[[1]]
h <- grep("^SAS name\\s+Counts", L) # page header ends here
if (length(h)) L <- L[-seq_len(h[1] + 1)]
L
}), use.names = FALSE)
}
nh3_codebook <- function(pdf, vars, lines = nh3_codebook_lines(pdf)) {
blank <- !nzchar(trimws(lines))
is_val <- function(x) grepl("^\\s*[0-9][0-9,]*\\s+\\S", x)
out <- lapply(vars, function(v) {
i <- grep(paste0("^\\s{0,8}", v, "(\\s|$)"), lines)[1]
if (is.na(i)) return(NULL)
## The entry ends at the first blank line after the value rows start
## (earlier blank lines separate the question from interviewer notes),
## or just before the next entry's position + name lines.
end <- i; seen <- FALSE; k <- i + 1
while (k <= length(lines)) {
if (grepl("^\\s{0,8}[A-Z][A-Z0-9_]*(\\s|$)", lines[k]) && k > i + 1 &&
grepl("^\\s*[0-9]+(-[0-9]+)?\\s", lines[k - 1])) { end <- k - 2; break }
if (blank[k]) { if (seen) break } else { end <- k; if (is_val(lines[k])) seen <- TRUE }
k <- k + 1
}
body <- if (end > i) lines[(i + 1):end] else character(0)
first <- which(is_val(body))[1]
pre <- if (!is.na(first) && first > 1) trimws(body[seq_len(first - 1)]) else character(0)
gap <- which(!nzchar(pre))[1]
if (is.na(gap)) gap <- length(pre) + 1
desc <- c(sub("^\\s*[0-9]+(-[0-9]+)?\\s*", "", lines[i - 1]),
sub(paste0("^\\s*", v, "\\s*"), "", lines[i]), pre[seq_len(gap - 1)])
desc <- sub("\\s{2,}See note\\s*$", "", trimws(desc))
notes <- pre[seq_along(pre) >= gap]
## the name line can itself carry the first value row (" HAG17A 418 1 Yes")
first_row <- sub(paste0("^\\s*", v), "", lines[i])
if (is_val(first_row)) { body <- c(first_row, body); desc <- desc[-2] }
m <- regmatches(body[is_val(body)],
regexec("^\\s*([0-9][0-9,]*)\\s+(\\S+)\\s*(.*?)\\s*$", body[is_val(body)], perl = TRUE))
tab <- data.frame(count = as.integer(gsub(",", "", vapply(m, `[`, "", 2))),
code = vapply(m, `[`, "", 3),
label = vapply(m, `[`, "", 4))
b <- tab$code == "Blank"; tab$label[b] <- "Blank"; tab$code[b] <- ""
list(var = v,
position = sub("^\\s*([0-9]+(-[0-9]+)?).*$", "\\1", lines[i - 1]),
description = paste(desc[nzchar(desc)], collapse = " "),
notes = paste(notes[nzchar(notes)], collapse = " "),
values = tab)
})
names(out) <- vars
if (any(lost <- vapply(out, is.null, logical(1))))
warning("not found: ", paste(vars[lost], collapse = ", "))
out[!lost]
}
## Codes made of a repeated 6, 7, 8 or 9 (666, 77, 8888 ...) -- each needs a
## decision before the variable is used as a number
nh3_special_codes <- function(cb) {
out <- do.call(rbind, lapply(cb, function(e) {
t <- e$values[grepl("^(6+|7+|8+|9+)$", e$values$code), ]
if (nrow(t)) data.frame(var = e$var, t, row.names = NULL)
}))
rownames(out) <- NULL
out
}
cb <- nh3_codebook(file.path(raw, "adult-acc.pdf"),
c("HAR1", "HAR3", "HAR4S", "HAR16", "HAR24", "HAR27"))
cb$HAR4S$values
#> count code label
#> 1 4656 001-140
#> 2 59 666 Varies or varied
#> 3 237 777 Less than 1 cigarette per day
#> 4 51 888 Blank but applicable
#> 5 5 999 Don't know
#> 6 15042 Blank
nh3_special_codes(cb)
#> var count code label
#> 1 HAR1 16 8 Blank but applicable
#> 2 HAR3 18 8 Blank but applicable
#> 3 HAR4S 59 666 Varies or varied
#> 4 HAR4S 237 777 Less than 1 cigarette per day
#> 5 HAR4S 51 888 Blank but applicable
#> 6 HAR4S 5 999 Don't know
#> ...For yes/no items such as HAR1 the table confirms the
coding (1 = yes, 2 = no, 8 = blank but applicable). A wrapped question
line that starts with a number (HAR23: “…at least / 20
cigars in your entire life?”) is misread as a value row, so glance at
the full entry for anything unexpected.
6. Smoking status, checked against cotinine
HAR1 asks whether the person has smoked 100 or more
cigarettes; only those who say yes are asked HAR3 (smoke
now). Other tobacco items are asked only of people who ever used that
product, so a skipped item means “no”.
d <- merge(exam, adult, by = "SEQN") # adults 17+ with an exam
d <- merge(d, lab2, by = "SEQN", all.x = TRUE)
d$cig_status <- with(d, ifelse(HAR1 == 2, "Never",
ifelse(HAR1 == 1 & HAR3 == 2, "Former",
ifelse(HAR1 == 1 & HAR3 == 1, "Current", NA))))
d$other_tob <- with(d, HAR16 %in% 1 | HAR24 %in% 1 | HAR27 %in% 1)
## serum cotinine > 10 ng/mL is the usual cutoff for active tobacco use
aggregate(COP ~ cig_status, data = d,
FUN = function(v) round(c(median = median(v), pct_over10 = 100 * mean(v > 10)), 2))About 6% of self-reported never-smokers have cotinine above 10 ng/mL, and about 13% of former smokers. A stricter never-smoker definition:
7. Link mortality
nhanesR downloads and parses the 2019 public-use Linked
Mortality File for NHANES III like any other cycle.
library(nhanesR)
mort <- nhanes_mortality_parse(cycles = "1988-1994")[["1988-1994"]]
d <- merge(d, mort[, c("SEQN", "ELIGSTAT", "MORTSTAT", "UCOD_LEADING",
"PERMTH_EXM", "PERMTH_INT")],
by = "SEQN", all.x = TRUE)
table(d$ELIGSTAT, d$MORTSTAT, useNA = "ifany")
d <- d[d$ELIGSTAT %in% 1 & !is.na(d$PERMTH_EXM), ]
d$time <- d$PERMTH_EXM # months from the exam
d$event <- d$MORTSTAT # 1 = died by 2019-12-31ELIGSTAT is 2 for people under 18 at the survey (not
linked) and 3 when there was too little identifying information. About
450 linkable adults have no PERMTH_EXM: they were examined
at home rather than in the mobile exam centre, and have a zero MEC exam
weight. PERMTH_INT gives their follow-up from the interview
if you need them.
8. Survey-weighted Cox model
library(survey)
library(survival)
d$Ht_m <- d$BMPHT / 100
des <- svydesign(ids = ~SDPPSU6, strata = ~SDPSTRA6, weights = ~WTPFEX6,
nest = TRUE, data = d[d$WTPFEX6 > 0, ])
fit <- svycoxph(Surv(time, event) ~ HSAGEIR + factor(HSSEX) + BMPBMI + Ht_m,
design = subset(des, never_conf & !is.na(BMPBMI)))
summary(fit)Notes
- Public-use mortality files add small perturbations to some follow-up times and causes of death for confidentiality. NCHS reports the effect on relative-risk estimates is negligible.
-
UCOD_LEADINGuses the same 10-category leading-cause recode as the continuous-NHANES files; seenhanes_ucod_labels(). - NHANES III also has youth (
1a/youth.dat) and first-laboratory (1a/lab.dat) files; the same reader works for them.