brfssdata is an R package, but the files it distributes are not an R format. Every survey year is a plain parquet file on a public release, and the R functions are one client among several. This article shows how to read those files from Python, from the DuckDB command line, and from SAS and Stata, and how to export an extract from R when a collaborator needs a file they can open in their own software.
The hosted files
Forty survey years, 1985 through 2024, are published as one parquet file per year at a fixed URL. Substitute the survey year in both places:
https://github.com/muntasirmasum/brfssdata/releases/download/data-YYYY/brfss_YYYY.parquetFor 2023 that resolves to:
https://github.com/muntasirmasum/brfssdata/releases/download/data-2023/brfss_2023.parquetThere is no token and no API. The 2023 file is about 29 MB and holds
433,323 rows and 351 columns, compressed with zstd, which every parquet
reader of the last few years handles. Columns keep CDC’s canonical
uppercase names, leading underscores on calculated variables included
(_LLCPWT, _STSTR, _PSU,
_IMPRACE), so a codebook written against CDC’s public-use
file applies unchanged. One column is added, year, which
matters once you read several files at once.
Almost every survey variable is stored as a physical DOUBLE,
mirroring the SAS numerics these files come from. The exceptions are
stored as strings, and which columns those are varies by year: the date
parts (IDATE, IMONTH, IDAY,
IYEAR) in every year, SEQNO from 1999 on,
INTVID through 2012, and a few text-coded columns in the
middle years (MRACE and RCSRACE among them).
The added year column is an integer. Storing the numerics
uniformly is deliberate and keeps a code comparable across files, but it
is not a guarantee about CDC’s own choices: a column published as text
in one year and as a number in another still differs between files,
which is what the brfssdata_type_conflict error reports
when such a pair is read together. The cost of DOUBLE storage is
cosmetic, since an integral code renders as 2.0 in a CSV
extract, so cast in SQL (GENHLTH::INTEGER) when integers
are wanted.
The metadata lives under a single data-meta tag, nine
assets in all. manifest.json lists the published years and
carries a per-asset sha256 and size. Four parquet catalogs sit beside
it: brfss_variables.parquet is the variable catalog (2,128
distinct variables with their labels and the years each one appears in),
brfss_labels.parquet is the value-label catalog covering
1998 onward, brfss_crosswalk.parquet holds CDC’s rename
families, so _DRNKWK1 through _DRNKWK3 read as
one concept with a review status and a comparability note, and
brfss_year_info.parquet is the per-year inventory of
respondents, variables, jurisdictions, file size, and codebook URL. Each
of the four catalogs has a .sha256 companion at the same
URL with .sha256 appended; manifest.json does
not, since it is the record the hashes live in. An air-gapped mirror
wants all nine. A mirror missing the crosswalk falls back to the
snapshot bundled with the installed package, which is only as current as
that package version, and a mirror missing the year inventory has no
bundled fallback to reach for at all.
https://github.com/muntasirmasum/brfssdata/releases/download/data-meta/manifest.json
https://github.com/muntasirmasum/brfssdata/releases/download/data-meta/brfss_variables.parquet
https://github.com/muntasirmasum/brfssdata/releases/download/data-meta/brfss_labels.parquet
https://github.com/muntasirmasum/brfssdata/releases/download/data-meta/brfss_crosswalk.parquet
https://github.com/muntasirmasum/brfssdata/releases/download/data-meta/brfss_year_info.parquetPython
pandas reads a parquet URL directly. Naming the columns you want keeps memory down, though on plain HTTP the whole file still arrives before anything is decoded.
import pandas as pd
url = (
"https://github.com/muntasirmasum/brfssdata/releases/download/"
"data-2023/brfss_2023.parquet"
)
df = pd.read_parquet(url, columns=["year", "GENHLTH", "_LLCPWT"])
df["GENHLTH"].value_counts().sort_index()DuckDB is the better tool when you want a slice of a year rather than
the year. It reads the file over HTTP range requests, so a query
touching three columns transfers roughly those three columns instead of
all 351. Quote the underscore-prefixed names in SQL. DuckDB itself
accepts them bare, since its lexer follows PostgreSQL in allowing a
leading underscore, so sum(_LLCPWT) runs; other engines do
not, and quoting is what makes the query portable.
import duckdb
duckdb.sql(f"""
SELECT GENHLTH,
count(*) AS n,
sum("_LLCPWT") AS weighted_n
FROM read_parquet('{url}')
GROUP BY GENHLTH
ORDER BY GENHLTH
""")Add .df() to hand the result back to pandas, or
.arrow() for an Arrow table. Several years combine in one
call. Variable sets drift across years, so
union_by_name = true fills the gaps with nulls, and the
year column tells the rows apart afterwards.
base = "https://github.com/muntasirmasum/brfssdata/releases/download"
files = [f"{base}/data-{y}/brfss_{y}.parquet" for y in (2021, 2022, 2023)]
duckdb.sql(f"""
SELECT year, count(*) AS n
FROM read_parquet({files}, union_by_name = true)
GROUP BY year
ORDER BY year
""").df()polars works the same way through pl.read_parquet(), and
pyarrow through pyarrow.parquet.read_table().
The DuckDB command line
Nothing about the above needs Python. The DuckDB CLI runs the same
queries, and the httpfs extension is what teaches it to
read an https path. Recent versions load that extension on demand, but
installing it explicitly works on any version.
INSTALL httpfs;
LOAD httpfs;
SELECT year, GENHLTH, count(*) AS n
FROM read_parquet(
'https://github.com/muntasirmasum/brfssdata/releases/download/data-2023/brfss_2023.parquet'
)
GROUP BY year, GENHLTH
ORDER BY GENHLTH;DESCRIBE SELECT * FROM read_parquet(...) prints the
schema without reading the rows, which is a quick way to check whether a
variable exists in a given year. To hand a subset to software that
cannot read parquet, COPY writes it out:
SAS
SAS support for parquet depends on which SAS you run, and the honest summary is narrow. The SAS Viya platform provides a LIBNAME engine for parquet, introduced in Viya 2021.2.6, and SAS documents ZSTD, the codec these files use, among the compression types that engine reads. That engine is not part of SAS 9, and its reach has grown release by release, starting with directories local to the compute server and adding cloud object storage later. It also points at a directory rather than a URL, so you would download the file first. If you run Viya, check the LIBNAME engines for ORC and Parquet reference for your release before assuming a path works; the shape of the code is this:
libname brfss parquet "/path/to/downloaded/files";
proc surveyfreq data = brfss.brfss_2023;
strata _STSTR;
cluster _PSU;
weight _LLCPWT;
tables GENHLTH;
run;The route that works on any SAS, including SAS 9, is to export a SAS
Transport file from R. haven::write_xpt() writes one.
library(brfssdata)
library(dplyr)
library(haven)
extract <- read_brfss(
2023,
vars = c(
"GENHLTH", "SEXVAR", "_AGE_G", "_IMPRACE",
"_LLCPWT", "_STSTR", "_PSU"
),
quiet = TRUE
) |>
slice_head(n = 5000)
xpt <- tempfile(fileext = ".xpt")
write_xpt(extract, xpt)
file.size(xpt)
#> [1] 321840Reading it back confirms the names survived:
names(read_xpt(xpt, n_max = 1))
#> [1] "GENHLTH" "SEXVAR" "_AGE_G" "_IMPRACE" "_LLCPWT" "_STSTR" "_PSU"
#> [8] "year"Which transport version
Transport files come in two flavors, SAS reads them by different
routes, and version is therefore the argument that decides
whether your collaborator can open the file at all.
Version 8 is haven’s default and allows variable names up to 32
characters. SAS’s XPORT libname engine cannot read it, because that
engine implements the Version 5 feature set; pointed at a V8 file it
reports that the file is not a SAS data set. What reads it is the
%XPT2LOC autocall macro, which has shipped with SAS since
9.4M2 and handles both flavors:
filename xptfile "/path/to/brfss23.xpt";
%xpt2loc(libref = work, filespec = xptfile);Version 5 is the older specification and the one the XPORT engine reads directly, which is why submission workflows and older sites still ask for it:
libname xptfile xport "/path/to/brfss23.xpt";
proc copy in = xptfile out = work;
run;Version 5 caps both the dataset member name and every variable name at 8 characters. The member name defaults to the file name without its extension, so a random temporary file name fails outright:
write_xpt(extract, tempfile(fileext = ".xpt"), version = 5)
#> Error in `write_xpt()`:
#> ! `name` must be 8 characters or fewer.Pass name to fix that. Variable names are the subtler
problem. Every one of the 2,128 variables in CDC’s catalog is 8
characters or fewer, _IMPRACE sitting exactly at the limit,
so a straight CDC extract comes through version 5 unchanged. Derived
variables do not, and the truncation is silent:
derived <- extract |>
mutate(fairpoor_health = as.integer(GENHLTH %in% 4:5))
v5 <- tempfile(fileext = ".xpt")
write_xpt(derived, v5, version = 5, name = "brfss23")
names(read_xpt(v5, n_max = 1))
#> [1] "GENHLTH" "SEXVAR" "_AGE_G" "_IMPRACE" "_LLCPWT" "_STSTR" "_PSU"
#> [8] "year" "fairpoor"fairpoor_health arrives as fairpoor, with
no warning, and two derived names sharing their first 8 characters would
collide the same way. Choose by what the recipient will run: version 8
if they can call %XPT2LOC, version 5 if they will point the
XPORT engine at the file, and in that second case rename your derived
variables to 8 characters yourself so you pick the abbreviations instead
of inheriting them.
Transport files carry the numeric codes, not the meanings. Ship the relevant slice of the value-label catalog next to the data so the recipient has a codebook:
brfss_labels(
c("GENHLTH", "SEXVAR", "_AGE_G", "_IMPRACE"),
years = 2023
) |>
arrange(variable, code)
#> # A tibble: 21 × 5
#> year variable code label complete
#> <int> <chr> <int> <chr> <lgl>
#> 1 2023 GENHLTH 1 Excellent TRUE
#> 2 2023 GENHLTH 2 Very good TRUE
#> 3 2023 GENHLTH 3 Good TRUE
#> 4 2023 GENHLTH 4 Fair TRUE
#> 5 2023 GENHLTH 5 Poor TRUE
#> 6 2023 GENHLTH 7 Dont know/Not Sure TRUE
#> 7 2023 GENHLTH 9 Refused TRUE
#> 8 2023 SEXVAR 1 Male TRUE
#> 9 2023 SEXVAR 2 Female TRUE
#> 10 2023 _AGE_G 1 Age 18 to 24 TRUE
#> # ℹ 11 more rowsStata
haven::write_dta() writes a Stata file, defaulting to
the Stata 14 format; pass version = 13 or lower for older
Stata. The complication is naming. Stata permits a leading underscore
but reserves such names for its own system variables (_n,
_N, _b, _cons), and the manual
advises against user variables that begin with one. haven does not
enforce this, so write_dta() on a raw BRFSS extract
succeeds and the trouble surfaces later, in Stata. Rename first:
for_stata <- extract |>
rename_with(\(x) sub("^_", "x_", x))
names(for_stata)
#> [1] "GENHLTH" "SEXVAR" "x_AGE_G" "x_IMPRACE" "x_LLCPWT" "x_STSTR"
#> [7] "x_PSU" "year"
dta <- tempfile(fileext = ".dta")
write_dta(for_stata, dta)
file.size(dta)
#> [1] 305675Prefixing is safer than deleting the underscore outright. CDC’s
catalog contains both _RACE and RACE, and both
_SEX and SEX, among other pairs, so stripping
the character can quietly merge two different variables in a multi-year
extract. Stata allows names up to 32 characters, which nothing in BRFSS
comes near.
On the Stata side, declare the design before estimating anything:
use "brfss_2023_extract.dta", clear
svyset x_PSU [pweight = x_LLCPWT], strata(x_STSTR) singleunit(centered)
svy: proportion GENHLTHSingle-PSU strata are common in BRFSS, because the public-use files
give each respondent their own primary sampling unit within a stratum,
so any stratum holding one respondent holds one PSU. Software differs in
what it does with them. In R, brfss_design() sets
options(survey.lonely.psu = "adjust") when the design it
just built actually contains such a stratum and the option is unset (any
value other than survey’s own load-time "fail" counts as
set by you), which centers those strata at the grand mean;
singleunit(centered) is the Stata option with the same
behavior. A year with no single-PSU stratum leaves the option alone, so
an unrelated survey analysis later in the session keeps survey’s own
default.
Value labels without losing CDC’s codes
read_brfss(labels = TRUE) returns factors for variables
whose CDC format is a complete code-to-label map, and
write_dta() does turn those factors into Stata value
labels. What it cannot turn them into is CDC’s codes. A factor carries
its levels and their positions, nothing else, so the export numbers the
levels 1, 2, 3 in order and the original codes are gone:
labeled <- read_brfss(2023, vars = "GENHLTH", labels = TRUE, quiet = TRUE)
factor_dta <- tempfile(fileext = ".dta")
write_dta(labeled, factor_dta)
read_dta(factor_dta) |>
count(value = as.numeric(GENHLTH), label = as_factor(GENHLTH))
#> # A tibble: 8 × 3
#> value label n
#> <dbl> <fct> <int>
#> 1 1 Excellent 63410
#> 2 2 Very good 142115
#> 3 3 Good 144209
#> 4 4 Fair 61955
#> 5 5 Poor 20372
#> 6 6 Dont know/Not Sure 897
#> 7 7 Refused 361
#> 8 NA NA 4“Dont know/Not Sure” arrives on value 6 and “Refused” on 7, where CDC
wrote 7 and 9. Reading with na = TRUE does not rescue this;
it drops the missing-type levels and renumbers whatever is left, which
moves substantive codes too. CHECKUP1’s code 8 is “Never”,
a real answer from 2,747 respondents in 2023, and under
na = TRUE it arrives on value 5.
The renumbering is silent, and Stata is where it does its damage,
because the reflex there is a line like
replace GENHLTH = . if inlist(GENHLTH, 7, 9). Against a
renumbered file that line matches nothing, and every don’t-know answer
stays in the estimation sample.
Two routes keep the codes. The plain one is to export the raw codes,
as for_stata does above, and ship the relevant slice of
brfss_labels() alongside as a codebook. The other builds
Stata value labels from that same catalog, so the .dta
carries the true codes and their meanings together:
codes <- brfss_labels(c("GENHLTH", "SEXVAR"), years = 2023)
label_from_catalog <- function(x, name) {
map <- codes[codes$variable == name, ]
map <- map[order(map$code), ]
labelled(x, setNames(map$code, map$label))
}
labeled_codes <- for_stata |>
mutate(
GENHLTH = label_from_catalog(GENHLTH, "GENHLTH"),
SEXVAR = label_from_catalog(SEXVAR, "SEXVAR")
)
codes_dta <- tempfile(fileext = ".dta")
write_dta(labeled_codes, codes_dta)
read_dta(codes_dta) |>
count(value = as.numeric(GENHLTH), label = as_factor(GENHLTH))
#> # A tibble: 7 × 3
#> value label n
#> <dbl> <fct> <int>
#> 1 1 Excellent 582
#> 2 2 Very good 1442
#> 3 3 Good 1752
#> 4 4 Fair 925
#> 5 5 Poor 286
#> 6 7 Dont know/Not Sure 8
#> 7 9 Refused 5CDC’s 7 and 9 survive, and label list in Stata shows the
wording against them. Do this only for the variables the analysis
touches. The catalog is per year and per variable, so a multi-year
extract needs its labels checked for drift before one year’s map is
applied to another.
CSV
CSV is the format nothing refuses, and when you do not know what the recipient will open the file with, it is the safe answer.
readr::write_csv(for_stata, "brfss_2023_extract.csv")
# base R only
utils::write.csv(for_stata, "brfss_2023_extract.csv", row.names = FALSE)The cost is size and types. Written out in full, the 2023 file is about 447 MB as CSV against 29 MB as parquet, roughly fifteen times larger. None of the column types survive the trip either, so every reader on the other end re-guesses whether a column is integer, double, or text, and a leading-zero FIPS code or a variable that is numeric in one year and character in another is where that guessing goes wrong. Send a column-selected, row-filtered extract instead of a whole year, and if the recipient can read parquet at all, send parquet.
Citing and verifying
Whatever tool opens the file, the data are CDC’s and the citation is CDC’s:
Centers for Disease Control and Prevention (CDC). Behavioral Risk Factor Surveillance System Survey Data. Atlanta, Georgia: U.S. Department of Health and Human Services, Centers for Disease Control and Prevention, [appropriate year].
Each parquet file is built from CDC’s published SAS Transport release
by the pipeline in data-raw/,
and every published survey year has a .sha256 companion at
the same URL with .sha256 appended. Anyone downloading
outside R can check what they received:
base=https://github.com/muntasirmasum/brfssdata/releases/download/data-2023
curl -LO $base/brfss_2023.parquet
curl -LO $base/brfss_2023.parquet.sha256
shasum -a 256 -c brfss_2023.parquet.sha256 # sha256sum -c on Linux
#> brfss_2023.parquet: OKRecording that checksum alongside the analysis is what lets someone else confirm years later that they are holding the same file you analyzed.