Using R with NHANES
Cubicle 3253
---
title: Self-Reported Snoring in Adults
subtitle: Using R with NHANES
author: Cubicle 3253
format:
html:
code-tools: true
echo: false
embed-resources: true
fig-asp: 0.6
fig-dpi: 300
grid:
sidebar-width: 200px
body-width: 1500px
margin-width: 100px
linestretch: 1.1
toc: true
toc-location: left
---
```{r}
#| message: false
library(apexcharter)
library(dplyr)
library(gt)
library(haven)
library(survey)
library(srvyr)
library(viridisLite)
read_nhanes <- function(xpt_name) {
xpt_name <- toupper(xpt_name)
if(substr(xpt_name, 1, 2) == "P_") {
begin_year = "2017"
} else {
begin_year = switch(substr(xpt_name, nchar(xpt_name) - 1, nchar(xpt_name)),
"_L" = "2021",
"_J" = "2017",
"_I" = "2015",
"_H" = "2013",
"_G" = "2011",
"_F" = "2009",
"_E" = "2007",
"_D" = "2005",
"_C" = "2003",
"_B" = "2001",
"1999")
}
xpt_url <- paste0("https://wwwn.cdc.gov/nchs/data/nhanes/public/", begin_year, "/datafiles/", xpt_name, ".xpt")
try(read_xpt(xpt_url), silent = TRUE)
}
```
```{r}
demo_vars <- c("SEQN", "SDDSRVYR", "RIAGENDR", "RIDAGEYR", "DMDMARTL", "DMDMARTZ",
"SDMVSTRA", "SDMVPSU", "WTINT2YR", "WTINTPRP")
DEMO_D <- read_nhanes("DEMO_D") |> select(any_of(demo_vars))
DEMO_E <- read_nhanes("DEMO_E") |> select(any_of(demo_vars))
DEMO_I <- read_nhanes("DEMO_I") |> select(any_of(demo_vars))
P_DEMO <- read_nhanes("P_DEMO") |> select(any_of(demo_vars))
DEMO <- bind_rows(DEMO_D, DEMO_E, DEMO_I, P_DEMO)
```
```{r}
SLQ_D <- read_nhanes("SLQ_D") |> select(SEQN, SLQ030)
SLQ_E <- read_nhanes("SLQ_E") |> select(SEQN, SLQ030)
SLQ_I <- read_nhanes("SLQ_I") |> select(SEQN, SLQ030)
P_SLQ <- read_nhanes("P_SLQ") |> select(SEQN, SLQ030)
SLQ <- bind_rows(SLQ_D, SLQ_E, SLQ_I, P_SLQ)
```
```{r}
One <- left_join(DEMO, SLQ, by = "SEQN") |>
mutate(
RIAGENDR = if_else(RIAGENDR == 1, "Men", "Women"),
RIDAGEYR = case_match(RIDAGEYR,
20:34 ~ "20 to 34",
35:49 ~ "35 to 49",
50:64 ~ "50 to 64",
65:80 ~ "65 or more"),
DMDMARITAL = case_when(
SDDSRVYR != 66 & DMDMARTL %in% c(1, 6) ~ "Married/Living with Partner",
SDDSRVYR != 66 & DMDMARTL %in% 2:4 ~ "Widowed/Divorced/Separated",
SDDSRVYR != 66 & DMDMARTL == 5 ~ "Never Married",
DMDMARTZ == 1 ~ "Married/Living with Partner",
DMDMARTZ == 2 ~ "Widowed/Divorced/Separated",
DMDMARTZ == 3 ~ "Never Married"),
DMDMARITAL = factor(DMDMARITAL,
levels = c("Married/Living with Partner", "Widowed/Divorced/Separated", "Never Married")),
## See: https://wwwn.cdc.gov/nchs/nhanes/continuousnhanes/overviewbrief.aspx?Cycle=2017-2020 ##
survey_wt = if_else(SDDSRVYR != 66, WTINT2YR * 2/9.2, WTINTPRP * 3.2/9.2),
SDDSRVYR = case_match(SDDSRVYR,
4 ~ "2005-2006",
5 ~ "2007-2008",
9 ~ "2015-2016",
66 ~ "2017-March 2020"),
SLQ030 = case_match(SLQ030,
0 ~ "Never",
1 ~ "Rarely",
2 ~ "Occasionally",
3 ~ "Frequently",
7:9 ~ "Refused/DK"),
SLQ030 = factor(SLQ030,
levels = c("Never", "Rarely", "Occasionally", "Frequently", "Refused/DK")))
```
```{r}
NHANES <- One |>
filter(!is.na(RIDAGEYR), !is.na(SLQ030)) |>
as_survey_design(id = SDMVPSU, strata = SDMVSTRA, nest = TRUE, weights = survey_wt)
```
## By Survey
<h4>Self-Reported Snoring in Adults by Survey Cycle</h4>
```{r}
t <- NHANES |>
group_by(SDDSRVYR, SLQ030) |>
summarize(p = survey_prop(proportion = TRUE, vartype = "ci", prop_method = "beta")) |>
mutate(across(p:p_upp, ~ .x * 100))
```
```{r}
#| results: asis
t |>
apex(aes(x = SDDSRVYR, y = p, fill = SLQ030), type = "bar", height = 540) |>
ax_chart(stacked = TRUE) |>
ax_colors(viridis(5, alpha = 0.8)) |>
ax_xaxis(title = list(text = "Percent", style = list(fontSize = "16px", fontWeight = 400)),
min = 0, max = 100, stepSize = 10,
labels = list(style = list(fontSize = "14px"))) |>
ax_yaxis(labels = list(align = "left", style = list(fontSize = "14px"))) |>
ax_plotOptions(bar = bar_opts(barHeight = "80%")) |>
ax_legend(position = "right") |>
ax_tooltip(y = list(formatter = format_num(".1f", suffix = "%")))
```
## By Sex & Survey
<h4>Self-Reported Snoring in Adults by Sex & Survey Cycle</h4>
```{r}
t <- NHANES |>
group_by(RIAGENDR, SDDSRVYR, SLQ030) |>
summarize(p = survey_prop(proportion = TRUE, vartype = "ci", prop_method = "beta")) |>
mutate(across(p:p_upp, ~ .x * 100))
```
```{r}
#| results: asis
t |>
apex(aes(x = SDDSRVYR, y = p, fill = SLQ030), type = "bar") |>
ax_chart(stacked = TRUE) |>
ax_colors(viridis(5, alpha = 0.8)) |>
ax_xaxis(title = list(text = "Percent", style = list(fontSize = "16px", fontWeight = 400)),
min = 0, max = 100, stepSize = 10,
labels = list(style = list(fontSize = "14px"))) |>
ax_yaxis(labels = list(align = "left", style = list(fontSize = "14px"))) |>
ax_plotOptions(bar = bar_opts(barHeight = "80%")) |>
ax_legend(position = "right") |>
ax_tooltip(y = list(formatter = format_num(".1f", suffix = "%"))) |>
ax_facet_wrap(vars(RIAGENDR), ncol = 1)
```
## By Age & Survey
<h4>Self-Reported Snoring in Adults by Age & Survey Cycle</h4>
```{r}
t <- NHANES |>
group_by(RIDAGEYR, SDDSRVYR, SLQ030) |>
summarize(p = survey_prop(proportion = TRUE, vartype = "ci", prop_method = "beta")) |>
mutate(across(p:p_upp, ~ .x * 100))
```
```{r}
#| results: asis
t |>
apex(aes(x = SDDSRVYR, y = p, fill = SLQ030), type = "bar", height = 150) |>
ax_chart(stacked = TRUE) |>
ax_colors(viridis(5, alpha = 0.8)) |>
ax_xaxis(title = list(text = "Percent", style = list(fontSize = "16px", fontWeight = 400)),
min = 0, max = 100, stepSize = 10,
labels = list(style = list(fontSize = "14px"))) |>
ax_yaxis(labels = list(align = "left", style = list(fontSize = "14px"))) |>
ax_plotOptions(bar = bar_opts(barHeight = "80%")) |>
ax_legend(position = "right") |>
ax_tooltip(y = list(formatter = format_num(".1f", suffix = "%"))) |>
ax_facet_wrap(vars(RIDAGEYR), ncol = 1)
```
## By Marital Status & Survey
<h4>Self-Reported Snoring in Adults by Age & Survey Cycle</h4>
```{r}
t <- NHANES |>
filter(!is.na(DMDMARITAL)) |>
group_by(DMDMARITAL, SDDSRVYR, SLQ030) |>
summarize(p = survey_prop(proportion = TRUE, vartype = "ci", prop_method = "beta")) |>
mutate(across(p:p_upp, ~ .x * 100))
```
```{r}
#| results: asis
t |>
apex(aes(x = SDDSRVYR, y = p, fill = SLQ030), type = "bar", height = 150) |>
ax_chart(stacked = TRUE) |>
ax_colors(viridis(5, alpha = 0.8)) |>
ax_xaxis(title = list(text = "Percent", style = list(fontSize = "16px", fontWeight = 400)),
min = 0, max = 100, stepSize = 10,
labels = list(style = list(fontSize = "14px"))) |>
ax_yaxis(labels = list(align = "left", style = list(fontSize = "14px"))) |>
ax_plotOptions(bar = bar_opts(barHeight = "80%")) |>
ax_legend(position = "right") |>
ax_tooltip(y = list(formatter = format_num(".1f", suffix = "%"))) |>
ax_facet_wrap(vars(DMDMARITAL), ncol = 1)
```