#' ------------------------------------------------------------------------------------------------------------
#' Example R code to replicate NCHS Data Brief No.303, Figures 1
#' Prevalence of Depression Among Adults Aged 20 and Over: United States, 2013-2016
#' Brody DJ, Pratt LA, Hughes JP. Prevalence of Depression Among Adults Aged 20 and Over: United
#' States, 2013-2016. NCHS Data Brief. No 303. Hyattsville, MD: National Center for Health Statistics. 2018.
#' Available at: https://www.cdc.gov/nchs/products/databriefs/db303.htm
#' ------------------------------------------------------------------------------------------------------------
# Load survey and dplyr packages
#+ message = FALSE, warning=FALSE
library(dplyr)
library(survey)
#'
options(survey.lonely.psu='adjust')
# Display Version Information
cat("R package versions:\n")
for (p in c("base", "survey","dplyr")) {
cat(p, ": ", as.character(packageVersion(p)), "\n")
}
#' # Data preparation
# Download & Read SAS Transport Files
# Demographic (DEMO)
download.file("https://wwwn.cdc.gov/Nchs/Data/Nhanes/Public/2013/DataFiles/DEMO_H.XPT", tf <- tempfile(), mode="wb")
DEMO_H <- foreign::read.xport(tf)[,c("SEQN","RIAGENDR","RIDAGEYR","SDMVSTRA","SDMVPSU","WTMEC2YR")]
download.file("https://wwwn.cdc.gov/Nchs/Data/Nhanes/Public/2015/DataFiles/DEMO_I.xpt", tf <- tempfile(), mode="wb")
DEMO_I <- foreign::read.xport(tf)[,c("SEQN","RIAGENDR","RIDAGEYR","SDMVSTRA","SDMVPSU","WTMEC2YR")]
# Mental Health - Depression Screener (DPQ)
download.file("https://wwwn.cdc.gov/Nchs/Data/Nhanes/Public/2013/DataFiles/DPQ_H.xpt", tf <- tempfile(), mode="wb")
DPQ_H <- foreign::read.xport(tf)
download.file("https://wwwn.cdc.gov/Nchs/Data/Nhanes/Public/2015/DataFiles/DPQ_I.xpt", tf <- tempfile(), mode="wb")
DPQ_I <- foreign::read.xport(tf)
# Append Files
DEMO <- bind_rows(DEMO_H, DEMO_I)
DPQ <- bind_rows(DPQ_H, DPQ_I)
# Merge DEMO and DPQ files and create derived variables
One <- left_join(DEMO, DPQ, by="SEQN") %>%
# Set 7=Refused and 9=Don't Know To Missing for variables DPQ010 thru DPQ090 ##
mutate_at(vars(DPQ010:DPQ090), ~ifelse(. >=7, NA, .)) %>%
mutate(. ,
# create indicator for overall summary
one = 1,
# Create depression score as sum of variables DPQ010 -- DPQ090
Depression.Score = rowSums(select(. , DPQ010:DPQ090)),
# Create depression indicator as binary 0/100 variable. (is missing if Depression.Score is missing)
Depression= ifelse(Depression.Score >=10, 100, 0),
# Create factor variables
Sex = factor(RIAGENDR, labels=c("Men", "Women")),
Age.Group = cut(RIDAGEYR, breaks=c(-Inf,19,39,59,Inf),labels=c("Under 20", "20-39","40-59","60 and over")),
# Generate 4-year MEC weight (Divide weight by 2 because we are appending 2 survey cycles)
# Note: using the MEC Exam Weights (WTMEC2YR), per the analytic notes on the
# Mental Health - Depression Screener (DPQ_H) documentation
WTMEC4YR = WTMEC2YR/2 ,
# Define indicator for analysis population of interest: adults aged 20 and over with a valid depression score
inAnalysis= (RIDAGEYR >= 20 & !is.na(Depression.Score))
) %>%
# drop DPQ variables
select(., -starts_with("DPQ"))
#' ## Define survey design
# Define survey design for overall dataset
NHANES_all <- svydesign(data=One, id=~SDMVPSU, strata=~SDMVSTRA, weights=~WTMEC4YR, nest=TRUE)
# Create a survey design object for the subset of interest: adults aged 20 and over with a valid depression score
# Subsetting the original survey design object ensures we keep the design information about the number of clusters and strata
NHANES <- subset(NHANES_all, inAnalysis)
#' ## Analysis
# Define a function to call svymean and unweighted count
getSummary <- function(varformula, byformula, design){
# Get mean, stderr, and unweighted sample size
c <- svyby(varformula, byformula, design, unwtd.count )
p <- svyby(varformula, byformula, design, svymean )
outSum <- left_join(select(c,-se), p)
outSum
}
#' ### Calculate prevalence of depression overall, by sex, by age group, and by age and sex
#' Adults
getSummary(~Depression, ~one, NHANES)
#' By sex
getSummary(~Depression, ~Sex, NHANES)
#' By age
getSummary(~Depression, ~Age.Group, NHANES)
#' By sex and age
getSummary(~Depression, ~Sex + Age.Group, NHANES)
#' ### Compare Prevalence Between Men And Women
svyttest(Depression~Sex, NHANES)$p.value %>% as.numeric
svyttest(Depression~Sex, subset(NHANES, Age.Group=="20-39"))$p.value %>% as.numeric
svyttest(Depression~Sex, subset(NHANES, Age.Group=="40-59"))$p.value %>% as.numeric
svyttest(Depression~Sex, subset(NHANES, Age.Group=="60 and over"))$p.value %>% as.numeric
#' ### Pairwise t-testing by age groups for total, men, and women
#' Differences by age group, among all adults
# Full output from svyttest command
svyttest(Depression~Age.Group, subset(NHANES, Age.Group=="20-39" | Age.Group=="40-59"))
# Displaying p-values only
svyttest(Depression~Age.Group, subset(NHANES, Age.Group=="20-39" | Age.Group=="40-59"))$p.value %>% as.numeric
svyttest(Depression~Age.Group, subset(NHANES, Age.Group=="20-39" | Age.Group=="60 and over"))$p.value %>% as.numeric
svyttest(Depression~Age.Group, subset(NHANES, Age.Group=="40-59" | Age.Group=="60 and over"))$p.value %>% as.numeric
#' Differences by age group, among men
svyttest(Depression~Age.Group, subset(NHANES, Sex=="Men" & (Age.Group=="20-39" | Age.Group=="40-59")))$p.value %>% as.numeric
svyttest(Depression~Age.Group, subset(NHANES, Sex=="Men" & (Age.Group=="20-39" | Age.Group=="60 and over")))$p.value %>% as.numeric
svyttest(Depression~Age.Group, subset(NHANES, Sex=="Men" & (Age.Group=="40-59" | Age.Group=="60 and over")))$p.value %>% as.numeric
#' Differences by age group, among women
svyttest(Depression~Age.Group, subset(NHANES, Sex=="Women" & (Age.Group=="20-39" | Age.Group=="40-59")))$p.value %>% as.numeric
svyttest(Depression~Age.Group, subset(NHANES, Sex=="Women" & (Age.Group=="20-39" | Age.Group=="60 and over")))$p.value %>% as.numeric
svyttest(Depression~Age.Group, subset(NHANES, Sex=="Women" & (Age.Group=="40-59" | Age.Group=="60 and over")))$p.value %>% as.numeric