--- title: "Getting started with brfssdata" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Getting started with brfssdata} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", eval = identical(Sys.getenv("IN_PKGDOWN"), "true") ) ``` The Behavioral Risk Factor Surveillance System (BRFSS) is the CDC's state-based telephone health survey. CDC established it in 1984, collects data in all 50 states, the District of Columbia, and participating US territories, and completes more than 400,000 adult interviews each year, which makes it the largest continuously conducted health survey system in the world. The public-use microdata are released one file per survey year, each carrying several hundred variables on health status, chronic conditions, health-care access, and risk behaviors. Getting at those files has traditionally meant downloading a zipped SAS transport archive for each year, reading it with a format-specific importer, and reconciling variable names by hand across years. brfssdata removes that step. Forty survey years, 1985 through 2024, are stored as compact parquet files on public releases; a year is downloaded once into a local cache and read from there afterwards, and every query runs through DuckDB, so asking for two variables out of a three-hundred-column survey reads two columns instead of the whole table. The installed vignette leaves most of this code unevaluated so the package builds offline, which means the copy you are reading in R shows no output under the code. The same document runs live against real data at , and reading it there is the better experience. ## Installation Once the package is on CRAN, install it the usual way. ```{r install, eval = FALSE} install.packages("brfssdata") ``` The development version comes from GitHub. ```{r install-dev, eval = FALSE} pak::pak("muntasirmasum/brfssdata") ``` ## A first read `read_brfss()` is the main entry point. Give it a survey year and the variables you want, and it returns a tibble with one row per respondent. ```{r first-read, message = FALSE} library(brfssdata) library(dplyr) dat <- read_brfss(2023, vars = c("GENHLTH", "PHYSHLTH")) dat ``` `GENHLTH` is self-rated general health on a five-point scale, and `PHYSHLTH` is the number of days in the past thirty on which physical health was not good. Both arrive as numeric codes, which is what most BRFSS analyses work from. A quick tabulation confirms the scale. ```{r first-count} dat |> count(GENHLTH) |> arrange(GENHLTH) ``` Leaving `vars` at its default (`NULL`) returns every column in the year, which for a recent release means well over three hundred variables. That is occasionally what you want, but naming the handful of variables you actually need is much faster and keeps memory use modest. ## Variable names BRFSS variable names are short, uppercase, and set by CDC. Names that begin with an underscore, such as `_LLCPWT`, `_STSTR`, or `_AGE_G`, are CDC's calculated variables: sampling weights, design identifiers, recoded age and race groupings, and derived risk-factor indicators. Everything else comes straight off the questionnaire. brfssdata always returns CDC's canonical spelling. An underscore-prefixed name is not a syntactic R name, so it needs backticks in dplyr code or brackets in base R (`` `_LLCPWT` `` or `dat[["_LLCPWT"]]`). Keeping CDC's spelling is deliberate. `_EDUCAG` is the name printed in the codebook, in CDC's own SAS and Stata examples, and in most published BRFSS code, so a column that still carries it can be looked up directly. Backticks are accepted everywhere a column name can appear, model formulas included, so a model fits straight on CDC's names. ```{r formula-backticks, eval = FALSE} svyglm(fairpoor ~ `_EDUCAG`, design = des, family = quasibinomial()) ``` If you would rather work in syntactic names throughout, `janitor::clean_names()` converts a whole tibble in one step. ```{r janitor, eval = FALSE} read_brfss(2023, vars = c("GENHLTH", "_EDUCAG")) |> janitor::clean_names() ``` The underscore is dropped rather than escaped, so `_EDUCAG` becomes `educag` and no longer matches the name you would search for in the codebook. There is a subtler cost as well. CDC ships both an underscored and a plain version of eight variables, and three of those pairs occupy the same survey years: `MRACE` and `_MRACE` both run from 2001 through 2012, `CHOLCHK` and `_CHOLCHK` overlap across many years between 1991 and 2015, and `FLUSHOT` and `_FLUSHOT` share 2002 and 2003. Cleaning the names makes each pair collide, and `clean_names()` settles the tie by suffixing whichever column happens to come second, so the two arrive as `mrace` and `mrace_2` with the assignment following column order rather than meaning. Check which is which before using either. Typing uppercase is tedious, so `vars` is matched case-insensitively. The lowercase request below finds `GENHLTH` and the returned column carries the canonical name. ```{r case-insensitive} names(read_brfss(2023, vars = "genhlth")) ``` Notice the `year` column. It is added to every result, whether or not you ask for it, so a tibble always knows which survey year each row came from. A name that matches nothing raises an error rather than silently returning fewer columns than you asked for. ```{r unknown-var, error = TRUE} read_brfss(2023, vars = "SMOKING") ``` ## Which years are available `brfss_years()` reads the manifest that accompanies the data releases and reports the published survey years. The manifest is cached and refreshed at most once a day, and `refresh = TRUE` forces a fresh copy. ```{r years} brfss_years() ``` ## Combining years Pass a vector of years and they come back stacked, with the `year` column identifying each. The complication is that BRFSS variable sets drift. Questions rotate on and off the core questionnaire, optional modules move between states, and CDC renames a calculated variable whenever its definition changes, usually by bumping a trailing digit. When a variable is absent from one of the requested years, it is filled with `NA` for those rows, and no error is raised. Computed weekly alcohol consumption is a clean example. The variable is `_DRNKWK1` in 2019 to 2021, `_DRNKWK2` in 2022 and 2023, and `_DRNKWK3` in 2024, with a change in definition at each rename. Asking for the first two across 2021 to 2023 shows where one gives way to the other. ```{r combine} drinks <- read_brfss(2021:2023, vars = c("_DRNKWK1", "_DRNKWK2")) drinks |> summarise( n = n(), drnkwk1 = sum(!is.na(`_DRNKWK1`)), drnkwk2 = sum(!is.na(`_DRNKWK2`)), .by = year ) |> arrange(year) ``` Reading both versions and coalescing them is the usual fix, but only after checking that the definitions are close enough to justify it. The package now tells you when you are standing in this trap: a request whose variable is empty in years a sibling generation covers gets a `brfssdata_rename_note` message pointing at `brfss_crosswalk()`, which lists the whole family (`_DRNKWK1` through `_DRNKWK3` here) along with its review status and, for reviewed pairs, whether consecutive generations stayed comparable and what changed. The decision to coalesce remains yours; the package hands you both columns and the family, and stays out of the recode. ## Finding variables Since names move, you generally want to search before you read. `brfss_vars()` matches a regular expression, case-insensitively, against both variable names and their labels, and reports the years each match appears in. ```{r search} brfss_vars("smok") ``` The labels being searched are CDC's SAS variable labels, capped at 40 characters and worded for the codebook rather than taken off the questionnaire. `PERSDOC3` carries `HAVE PERSONAL HEALTH CARE PROVIDER?`, so searching `doctor` in 2023 returns four other variables and not that one, and 135 of that year's labels sit at the cap with the last word cut in half. Search single words instead of phrases, try a synonym or two, and try the stem CDC would have used in a name, since names abbreviate (`doctor` reaches `PERSDOC3` only as `doc`). A search that finds nothing is more often a vocabulary mismatch than an absent variable. The `years` column collapses runs of consecutive years, so `2005-2024` means every year in that span and `2011-2013, 2020` means a variable that appeared for three years, disappeared, and came back once. Reading that column carefully is the fastest way to spot renamed variables. Two entries with near-identical labels and non-overlapping year ranges are almost always the same question under a new name. Calling `brfss_vars()` with no pattern returns the whole catalog, which holds 2,128 distinct variables across the forty years. Restricting to particular years narrows the search, and the result here shows the rename pattern plainly. ```{r search-years} brfss_vars("binge", years = 2021:2023) ``` `_RFBING5` and `_RFBING6` share a label and cover 2021 and then 2022 to 2023, which is the signature of a calculated variable whose definition was revised. The underlying questionnaire item, `DRNK3GE5`, kept its name across all three years. ## Value labels BRFSS answers are numeric codes, and the meaning of each code lives in CDC's SAS format libraries. `brfss_labels()` exposes those as a plain table with one row per year, variable, and code. ```{r labels-genhlth} brfss_labels("GENHLTH", years = 2023) ``` Setting `labels = TRUE` on a read converts eligible variables to factors using those maps. ```{r labels-true} read_brfss(2023, vars = c("GENHLTH", "SEXVAR"), labels = TRUE) ``` Two sex variables exist and they are not the same column. `SEXVAR` is the questionnaire response; `_SEX` is CDC's calculated version, the one the weighting is built on, and the two disagree for 178 of 2023's 433,323 records. Choose between them deliberately and name the choice in the write-up. Labeling renames codes; deciding which of them mean missing is a separate step, and the `GENHLTH` above has seven levels rather than five because CDC's don't-know and refused codes became factor levels like any other. The `na` argument handles that: `na = TRUE` sets the codes CDC uses for missing-type answers to `NA` before labeling, using the same catalog, so the factor arrives with its five substantive levels. ```{r labels-na} read_brfss(2023, vars = "GENHLTH", labels = TRUE, na = TRUE) |> count(GENHLTH) ``` `read_brfss()` defaults to `na = FALSE` and returns the file exactly as CDC published it; `brfss_design()` defaults to `na = TRUE`, because an estimate whose denominator includes "Refused" is almost never the estimate anyone wants. `brfss_missing_codes()` lists exactly which year, variable, and code combinations are affected, so the behavior is auditable rather than magic. ```{r missing-codes} brfss_missing_codes("GENHLTH", years = 2023) ``` Conversion is deliberately conservative. A variable becomes a factor only when its CDC format is a pure code-to-label map, carrying no numeric ranges such as `1-30` days, when the map is an unambiguous one-to-one correspondence (some CDC formats give one code several labels, or reuse one label across codes, and converting those would quietly rewrite the data), and when every value observed in the data falls inside that map. The `complete` column in the label catalog marks the formats that meet the first condition. `PHYSHLTH` is a good counterexample: its format documents only the special codes, because the substantive values are a count of days. ```{r labels-physhlth} brfss_labels("PHYSHLTH", years = 2023) ``` Turning that into a factor would silently destroy the day counts, so `labels = TRUE` leaves it numeric, and the catalog tells you that 88 means none, 77 means don't know, and 99 means refused. `na = TRUE` clears 77 and 99 here too, but 88 is an answer (zero days), not a missing code, and averaging `PHYSHLTH` without recoding it to 0 first overstates bad-health days for the healthiest respondents. ```{r physhlth-recode, eval = FALSE} read_brfss(2023, vars = "PHYSHLTH", na = TRUE) |> mutate(PHYSHLTH = replace(PHYSHLTH, PHYSHLTH == 88, 0)) ``` The same caution applies across years: when several years are requested, a variable converts only if its code set agrees across them. Labels cover 1998 onward. CDC does not distribute usable format libraries for earlier years, so `labels = TRUE` has nothing to work with before 1998 and quietly leaves those variables as numeric codes. ## Caching and offline work Every download lands in a per-user cache directory, resolved through `tools::R_user_dir()`. You can see where it is and what is in it. ```{r cache} brfss_cache_dir() brfss_cache_info() ``` Because reads come from the cache, a year you have already downloaded is available with no network at all. Setting `download = FALSE` makes that explicit: cached years are read normally, and a request for a year you do not have fails with a clear message instead of reaching for the network. ```{r offline} nrow(read_brfss(2023, vars = "GENHLTH", download = FALSE)) ``` This is the setting to use on a compute cluster without outbound network access, or in a reproducible pipeline where an unexpected download would be a bug. `brfss_download()` populates the cache in one call, years and metadata catalogs together, so the usual offline recipe is to run it once on a connected machine, copy the directory across, and point `options(brfssdata.cache_dir = ...)` at it. ```{r prefetch, eval = FALSE} brfss_download(2019:2023) ``` ### Restricted networks The data files live on GitHub releases, so any network that blocks or intercepts traffic to `github.com` blocks the download too. The common cases are hospital and corporate networks that terminate TLS at a proxy, campus networks with allowlists, and HPC compute nodes with no outbound access. A blocked download fails with an error naming the likely cause, not a hang, and nothing partial is written to the cache. The prefetch recipe above is the way through. Populate the cache from an unrestricted machine (a survey year from 2011 on is 20 to 45 MB, 2011 itself the largest, and all 40 years about 737 MB), copy the directory, and set `options(brfssdata.cache_dir = ...)` on the restricted one. A lab or a cluster needs only one such copy, since every user can point the same option at a shared directory. Variable discovery is never blocked. `brfss_vars()`, `brfss_labels()`, `brfss_codebook()`, and `brfss_crosswalk()` fall back to catalog snapshots bundled with the package, so codebook work runs before anything is downloaded. Downloads are verified against the sha256 checksums published with each release, and a cached file that later turns out damaged is either re-downloaded automatically or named in an error along with the exact call that fixes it. `brfss_cache_clear()` removes cached survey years, all of them when called with no arguments, only the ones you name otherwise. The manifest and catalogs stay unless you also pass `catalogs = TRUE`, so clearing data does not break offline variable searches. ```{r cache-clear, eval = FALSE} brfss_cache_clear(2021) brfss_cache_clear() brfss_cache_clear(catalogs = TRUE) ``` ## Survey-weighted analysis BRFSS uses a complex sampling design, and unweighted estimates from it are wrong in both the point estimate and the standard error. `brfss_design()` returns a [srvyr](https://cran.r-project.org/package=srvyr) design object with the sampling weight and strata (`_STSTR`) already applied, and the primary sampling units (`_PSU`) too in the years where those identify a real cluster, so you can move straight to analysis. The design states the specification it built, in `svyset` terms, as it is created. A design for 2001 or later prints `ids: 1`, and that is correct rather than a dropped design feature. CDC's public-use files from 2001 on number each respondent as their own PSU, so a cluster term would carry one level per row and change nothing: on 2023 the fair-or-poor `GENHLTH` estimate, its standard error, and its degrees of freedom are identical to the last bit with and without it. Through 2000 several respondents do share a `_PSU`, the clustering is real, and the design keeps it. ```{r design, message = FALSE} library(srvyr) brfss_design(2023, vars = "GENHLTH") |> filter(!is.na(GENHLTH)) |> group_by(GENHLTH) |> summarize(pct = survey_prop(vartype = "ci")) ``` The weight is chosen to match the survey era: `_FINALWT` for years before 2011 and `_LLCPWT` from 2011 on (a `weight` argument selects another final weight, such as the child weight `_CLLCPWT`, when CDC's documentation calls for it; a module weight like that exists only for its module's records, so the design subsets to the rows it covers and says how many were dropped; anything that is not a CDC final analysis weight, including intermediate pipeline stages like `_LLCPWT2`, is refused unless `unsafe_weight = TRUE` says you mean it; see `?brfss_design`). The design is built on three added columns, `brfss_wt`, `brfss_psu`, and `brfss_strata`, because CDC's names are not syntactic and cannot enter a model formula; the original CDC columns are kept alongside them, unchanged. The don't-know and refused codes arrive as `NA` here by default (`na = TRUE`), which is why the filter above is on `is.na()` and not on code ranges; pass `na = FALSE` for the raw codes. `brfss_design()` sets a few defaults that change results. Because BRFSS public-use files make each respondent their own primary sampling unit, single-PSU strata are common and would otherwise make variance estimation fail, so `brfss_design()` sets `options(survey.lonely.psu = "adjust")` if the option is unset, and says so once per session. Any value you set other than `"fail"` is respected as yours; `"fail"` is what the survey package itself installs on load, so it reads as unset, and `options(brfssdata.lonely_psu = ...)` is the way to pin any handling, `"fail"` included. When you request several years, weights are divided by the number of years, so pooled estimates describe an average year, not a summed population, while the strata become the year-by-stratum interaction, which treats each annual survey as an independent sample. Pass `pool_weights = FALSE` to leave the weights undivided; if state participation differs across the pooled years, a warning names the states involved. And when a requested variable has data almost only where a module weight such as `_CLLCPWT` does, a warning suggests that weight, because a module analysis under the full-sample default is very likely wrong; disable the check with `options(brfssdata.module_weight_check = FALSE)` if a state-optional module legitimately uses the core weight. ### The 2011 boundary In 2011 BRFSS combined cell phone and landline samples into a single public-use dataset for the first time, and adopted a new weighting method, iterative proportional fitting, also known as raking, in place of post-stratification. CDC advises data users not to make direct comparisons with data collected before 2011, and to begin new trend lines with that year. `brfss_design()` enforces the advice instead of burying it in a help page, so a request spanning the boundary fails. ```{r break-error, error = TRUE} brfss_design(2009:2013) ``` If you have a considered reason to pool across it, opt in explicitly and the call proceeds with a warning. ```{r break-allow, eval = FALSE} brfss_design(2009:2013, allow_break = TRUE) ``` Two articles carry this further. [Survey design in BRFSS](https://muntasirmasum.github.io/brfssdata/articles/survey-design.html) works through weights, strata, single-PSU handling, subpopulation analysis, and multi-year pooling in detail. [Survey-weighted logistic regression](https://muntasirmasum.github.io/brfssdata/articles/logistic-regression.html) fits a full model with `svyglm()`, walking from recoding and weighted prevalence through design-based confidence intervals on the odds ratios. ## Using the data outside R Nothing about the cache is specific to R. The files are ordinary parquet, readable by anything that speaks the format, so if you or a collaborator work in Python, SAS, Stata, or plain SQL, you can point those tools at the same files brfssdata downloaded, or export a prepared subset from R. [Using the data outside R](https://muntasirmasum.github.io/brfssdata/articles/outside-r.html) covers each route. ## Citing The package and the survey data are cited separately, and both entries come back from one call. ```{r citation, eval = TRUE} citation("brfssdata") ``` The second entry is CDC's recommended form for the data itself. Replace `[appropriate year]` with the survey year or years you analyzed, and repeat the citation for each year if your journal expects that. Manuscripts reporting BRFSS estimates conventionally also report response rates from CDC's Summary Data Quality Report for the years analyzed.