This document processes and analyzes data from major tobacco product use surveys that use complex survey design. The analysis involves calculating various survey metrics, including weighted counts, weighted percentages with 95% confidence intervals, and unweighted counts.
The sample data accompanying this sample code are for instructional purposes only. Please access the data from the original data providers for research or other purposes.
Steps include:
Each step is detailed with specific code snippets to achieve the desired analysis.
Population Assessment of Tobacco and Health (PATH) Study
# remove all the objects present in the workspace
rm(list=ls())
# load the data
# this sample data is a subset of PATH wave1
PATH_wave1 <- read_csv("sample_data/PATH_wave1_data.csv")
# These sample data are for instructional purposes only. This instructional use has been reviewed and approved by the ICPSR Data Stewardship Policy Committee.
# tobacco product use status
## cigarette
## never: never smoked 100+ cigarettes in their lifetime
## current: ever smoked 100+ cigarettes in their lifetime and used at least 1 day in the past 30 days
## former: neither "current" nor "never" users (complement)
## missing: missing data
table(PATH_wave1$UMp1_W1_A_smkstat_Reg_1past30, useNA = "always")##
## current former missing never <NA>
## 11208 5079 173 15849 0## e-cigarette
## never: never vaped e-cigarettes fairly regularly
## current: ever vaped e-cigarettes fairly regularly and used at least 1 day in the past 30 days
## former: neither "current" nor "never" users (complement)
## missing: missing data
table(PATH_wave1$UMp1_W1_A_ecigstat_Reg_1past30, useNA = "always")##
## current former missing never <NA>
## 1369 892 70 29978 0# age: 1: 18-24, 2: 25-34, 3: 35-54, 4: 55+
table(PATH_wave1$UM_W1_A_agecat, useNA = "always")##
## 1 2 3 4 <NA>
## 9110 6338 9778 7083 0# sex: 0: Female, 1: Male
table(PATH_wave1$UM_W1_A_male_imp, useNA = "always")##
## 0 1 <NA>
## 15989 16320 0# Filter out rows with missing value
PATH_wave1 <- PATH_wave1[!is.na(PATH_wave1$UM_W1_A_agecat),] # age
PATH_wave1 <- PATH_wave1[!is.na(PATH_wave1$UM_W1_A_male_imp),] # sex# generate a combination of tobacco product use status for cigarettes and e-cigarettes
# creates a vector that contains all possible combinations
cig_ecig_name<-paste0(rep(c("Never","Former","Current"),each=3),"/",
rep(c("Never","Former","Current")))
# loop through all combinations and assign values:
PATH_wave1$UM_W1_cig_ecig_9state<-NA
m<-1
for(i in c("never","former","current")){
for(j in c("never","former","current")){
PATH_wave1$UM_W1_cig_ecig_9state[PATH_wave1$UMp1_W1_A_smkstat_Reg_1past30 %in% i & PATH_wave1$UMp1_W1_A_ecigstat_Reg_1past30 %in% j]<-cig_ecig_name[m]
m<-m+1
}
}
# For example, the loop will generate and assign the combination "Current/Never" to the UM_W1_cig_ecig_9state column in PATH_wave1 for observations where cigarettes are "Current" and e-cigarettes are "Never".
PATH_wave1$UM_W1_cig_ecig_9state<-factor(PATH_wave1$UM_W1_cig_ecig_9state, levels=unique(cig_ecig_name))
table(PATH_wave1$UM_W1_cig_ecig_9state,useNA = "always")##
## Never/Never Never/Former Never/Current Former/Never Former/Former
## 15494 114 201 4536 187
## Former/Current Current/Never Current/Former Current/Current <NA>
## 351 9782 588 816 240PATH_wave1<-PATH_wave1[!is.na(PATH_wave1$UM_W1_cig_ecig_9state),]
table(PATH_wave1$UM_W1_cig_ecig_9state,useNA = "always")##
## Never/Never Never/Former Never/Current Former/Never Former/Former
## 15494 114 201 4536 187
## Former/Current Current/Never Current/Former Current/Current <NA>
## 351 9782 588 816 0# create a survey design object using replicate weights
svy_design <- svrepdesign(
id = ~PERSONID, # primary sampling unit (PSU) variable
weights = ~R01_A_PWGT, # sampling weights
repweights = "R01_A_PWGT[1-9]+", # replicate weights for PATH wave 1
type = "Fay", # replication method
rho = 0.3, # set the rho parameter for Fay's method
data = PATH_wave1 # your data set
)# calculate overall survey metrics by grouping the data by tobacco use status combinations
# compute weighted counts, percentages with 95% confidence intervals, and unweighted counts
# add overall indicators for sex and age
metrics_all <- as_survey_rep(svy_design) %>% # use as_survey_rep for replicate weights
# group by tobacco use status combinations, and ensures that all levels are included
group_by(UM_W1_cig_ecig_9state, .drop = FALSE) %>%
# calculate the outputs for each group
# please note that calculating the variance for the weighted percent can be quite slow for large complex data sets with many groups
summarize("weighted count" = survey_total(),
"weighted percent" = survey_mean(vartype = c("se","ci"), level = 0.95, proportion = TRUE),
"unweighted_count" = unweighted(n())) %>%
# indicate that the metrics are for which group in our later variables, here for the overall population
mutate(sex = "Overall",age = "Overall") %>%
ungroup()
metrics_all## # A tibble: 9 × 10
## UM_W1_cig_ecig_9state `weighted count` `weighted count_se` `weighted percent`
## <fct> <dbl> <dbl> <dbl>
## 1 Never/Never 143665118. 1257945. 0.613
## 2 Never/Former 356948. 43283. 0.00152
## 3 Never/Current 661093. 49619. 0.00282
## 4 Former/Never 45811205. 950333. 0.195
## 5 Former/Former 695678. 47417. 0.00297
## 6 Former/Current 1392067. 99235. 0.00594
## 7 Current/Never 37003647. 543252. 0.158
## 8 Current/Former 2017765. 105078. 0.00860
## 9 Current/Current 2915268. 131520. 0.0124
## # ℹ 6 more variables: `weighted percent_se` <dbl>,
## # `weighted percent_low` <dbl>, `weighted percent_upp` <dbl>,
## # unweighted_count <int>, sex <chr>, age <chr># calculation is similar to the overall survey metrics
# but it additionally groups the data by age categories, providing outputs for each age group within each tobacco use status combination
metrics_by_age <- as_survey_rep(svy_design) %>%
group_by(UM_W1_A_agecat, UM_W1_cig_ecig_9state, .drop = FALSE) %>%
summarize("weighted count" = survey_total(),
"weighted percent" = survey_mean(vartype = c("se", "ci"), level = 0.95, proportion = TRUE),
"unweighted_count" = unweighted(n())) %>%
mutate(sex = "Overall") %>%
ungroup() %>%
rename(age = UM_W1_A_agecat)
metrics_by_age## # A tibble: 36 × 10
## age UM_W1_cig_ecig_9state `weighted count` `weighted count_se`
## <dbl> <fct> <dbl> <dbl>
## 1 1 Never/Never 22714899. 213262.
## 2 1 Never/Former 175217. 24397.
## 3 1 Never/Current 276941. 32326.
## 4 1 Former/Never 1250580. 73082.
## 5 1 Former/Former 135716. 21147.
## 6 1 Former/Current 140173. 19486.
## 7 1 Current/Never 4875777. 156058.
## 8 1 Current/Former 441398. 31898.
## 9 1 Current/Current 549241. 45442.
## 10 2 Never/Never 25799126. 538644.
## # ℹ 26 more rows
## # ℹ 6 more variables: `weighted percent` <dbl>, `weighted percent_se` <dbl>,
## # `weighted percent_low` <dbl>, `weighted percent_upp` <dbl>,
## # unweighted_count <int>, sex <chr># calculation is similar to the overall survey metrics
# but it additionally groups the data by age and sex, providing outputs for each combination of sex and age within each tobacco use status combination
metrics_by_sex_age <- as_survey_rep(svy_design) %>%
group_by(UM_W1_A_male_imp, UM_W1_A_agecat, UM_W1_cig_ecig_9state, .drop = FALSE) %>%
summarize("weighted count" = survey_total(),
"weighted percent" = survey_mean(vartype = c("se", "ci"), level = 0.95, proportion = TRUE),
"unweighted_count" = unweighted(n())) %>%
ungroup() %>%
rename(sex = UM_W1_A_male_imp,age = UM_W1_A_agecat)
metrics_by_sex_age## # A tibble: 72 × 10
## sex age UM_W1_cig_ecig_9state `weighted count` `weighted count_se`
## <dbl> <dbl> <fct> <dbl> <dbl>
## 1 0 1 Never/Never 11924834. 127216.
## 2 0 1 Never/Former 67348. 12716.
## 3 0 1 Never/Current 99630. 17348.
## 4 0 1 Former/Never 627107. 52122.
## 5 0 1 Former/Former 60928. 14118.
## 6 0 1 Former/Current 50118. 11364.
## 7 0 1 Current/Never 2062679. 86122.
## 8 0 1 Current/Former 156809. 21357.
## 9 0 1 Current/Current 164508. 23048.
## 10 0 2 Never/Never 13784165. 350721.
## # ℹ 62 more rows
## # ℹ 5 more variables: `weighted percent` <dbl>, `weighted percent_se` <dbl>,
## # `weighted percent_low` <dbl>, `weighted percent_upp` <dbl>,
## # unweighted_count <int># combine survey metrics into a data frame
data <- rbind.data.frame(metrics_all, metrics_by_age, metrics_by_sex_age)
# select relevant columns and recode for readability
PATH_wave1_metrics <- data %>%
dplyr::select(UM_W1_cig_ecig_9state, age, sex,
"weighted count", "unweighted_count", "weighted percent",
"weighted percent_se", "weighted percent_low", "weighted percent_upp") %>%
# filters out rows with zero unweighted counts
filter(unweighted_count > 0)
# add wave identifier and rename columns
PATH_wave1_metrics <- PATH_wave1_metrics %>% mutate(wave = 1)
PATH_wave1_metrics <- rename(PATH_wave1_metrics, cig_cigar_ecig = UM_W1_cig_ecig_9state)
PATH_wave1_metrics <- rename(PATH_wave1_metrics, "unweighted count" = unweighted_count)
# recode sex and age variables
PATH_wave1_metrics$sex[PATH_wave1_metrics$sex %in% 0] <- "Female"
PATH_wave1_metrics$sex[PATH_wave1_metrics$sex %in% 1] <- "Male"
PATH_wave1_metrics$age[PATH_wave1_metrics$age %in% 1] <- "18-24"
PATH_wave1_metrics$age[PATH_wave1_metrics$age %in% 2] <- "25-34"
PATH_wave1_metrics$age[PATH_wave1_metrics$age %in% 3] <- "35-54"
PATH_wave1_metrics$age[PATH_wave1_metrics$age %in% 4] <- "55+"
# export the processed data
write.csv(PATH_wave1_metrics,"PATH_Wave1_metrics.csv",row.names = F)National Health Interview Survey (NHIS)
# remove all the objects present in the workspace
rm(list=ls())
# load the data
# this sample data is a subset of NHIS 2023
NHIS_2023 <- read_csv("sample_data/NHIS_2023_data.csv")
# These sample data are for instructional purposes only.
# tobacco product use status
## cigarette
## 1=current: ever smoked 100+ cigarettes in their lifetime and current every day/some day user
## 2=former: neither "current" nor "never" users (complement)
## 3=never: never smoked 100+ cigarettes in their lifetime
table(NHIS_2023$cig_status, useNA = "always")##
## 1 2 3 <NA>
## 3135 7235 18111 1041## e-cigarette
## 1=current: ever vaped e-cigarettes and current every day/some day user
## 2=former: neither "current" nor "never" users (complement)
## 3=never: never vaped e-cigarettes
table(NHIS_2023$ecig_status, useNA = "always")##
## 1 2 3 <NA>
## 1514 3518 23470 1020# Age: 1=18-24, 2=25-34, 3=35-44, 4=45-54, 5=55-64, 6=65+
table(NHIS_2023$agegrp, useNA = "always")##
## 1 2 3 4 5 6 <NA>
## 1857 4253 4540 4078 5053 9741 0# Sex: 0=Female, 1=Male
table(NHIS_2023$male, useNA = "always")##
## 0 1 <NA>
## 16059 13457 6# Filter out rows with missing value
NHIS_2023<-NHIS_2023[!is.na(NHIS_2023$agegrp),] # age
NHIS_2023<-NHIS_2023[!is.na(NHIS_2023$male),] # sex# generate a combination of tobacco product use status for cigarettes and e-cigarettes
# replace the numeric codes with the corresponding descriptive labels
NHIS_2023 <- NHIS_2023 %>%
mutate(
cig_status = recode(cig_status, `1` = "current", `2` = "former", `3` = "never"),
ecig_status = recode(ecig_status, `1` = "current", `2` = "former", `3` = "never")
)
# creates a vector that contains all possible combinations
cig_ecig_name<-paste0(rep(c("Never","Former","Current"),each=3),"/",c("Never","Former","Current"))
# loop through all combinations and assign values:
NHIS_2023$NHIS_9state<-NA
k<-1
for(i in c("never","former","current")){
for(j in c("never","former","current")){
NHIS_2023$NHIS_9state[NHIS_2023$cig_status %in% i & NHIS_2023$ecig_status %in% j]<-cig_ecig_name[k]
k<-k+1
}
}
# For example, the loop will generate and assign the combination "Current/Former" to the NHIS_9state column for observations where cigarettes are "Current" and e-cigarettes are "Former".
NHIS_2023$NHIS_9state<-factor(NHIS_2023$NHIS_9state, levels=unique(cig_ecig_name))
NHIS_2023<-NHIS_2023[!is.na(NHIS_2023$NHIS_9state),]
table(NHIS_2023$NHIS_9state,useNA = "always")##
## Never/Never Never/Former Never/Current Former/Never Former/Former
## 16415 1299 389 5311 1238
## Former/Current Current/Never Current/Former Current/Current <NA>
## 682 1718 975 439 0# set options for survey analysis
options(survey.lonely.psu = "adjust") # adjusts for lonely primary sampling units (PSUs) by redistributing weights
options(warn = -1) # suppress warnings
options(digits = 3) # shows 3 significant digits to display for numerical output
# create a survey design object using the svydesign function
svy_design <- svydesign(
id = ~fpx, # primary sampling unit (PSU) variable
strata = ~pstrat, # stratification variable
weights = ~wtfa_sa, # sampling weights
PSU = ~ppsu, # primary sampling unit (PSU)
data = NHIS_2023, # your data set
nest = TRUE # specify that the strata and PSUs are nested
)# calculate overall survey metrics by grouping the data by tobacco use status combinations
# compute weighted counts, percentages with 95% confidence intervals, and unweighted counts
# add overall indicators for sex and age
metrics_all <- as_survey(svy_design) %>% # convert the svy_design to a survey design object
# group by tobacco use status combinations, and ensures that all levels are included
group_by(NHIS_9state, .drop = FALSE) %>%
# calculate the outputs for each group
summarize("weighted count" = survey_total(),
"weighted percent" = survey_mean(vartype = c("se","ci"), level = 0.95, proportion = TRUE),
"unweighted_count" = unweighted(n())) %>%
# indicate that the metrics are for which group, here for the overall population
mutate(sex = "Overall",age = "Overall") %>%
ungroup()
metrics_all## # A tibble: 9 × 10
## NHIS_9state `weighted count` `weighted count_se` `weighted percent`
## <fct> <dbl> <dbl> <dbl>
## 1 Never/Never 146919837. 1045216. 0.591
## 2 Never/Former 14270696. 460667. 0.0574
## 3 Never/Current 4743082. 285991. 0.0191
## 4 Former/Never 37902940. 561263. 0.152
## 5 Former/Former 11011062. 357467. 0.0443
## 6 Former/Current 7000193. 311200. 0.0281
## 7 Current/Never 13837362. 386528. 0.0556
## 8 Current/Former 8585694. 320869. 0.0345
## 9 Current/Current 4511436. 248812. 0.0181
## # ℹ 6 more variables: `weighted percent_se` <dbl>,
## # `weighted percent_low` <dbl>, `weighted percent_upp` <dbl>,
## # unweighted_count <int>, sex <chr>, age <chr># calculation is similar to the overall survey metrics
# but it additionally groups the data by age categories, providing outputs for each age group within each tobacco use status combination
metrics_by_age <- as_survey(svy_design) %>%
group_by(agegrp, NHIS_9state, .drop = FALSE) %>%
summarize("weighted count" = survey_total(),
"weighted percent" = survey_mean(vartype = c("se", "ci"), level = 0.95, proportion = TRUE),
"unweighted_count" = unweighted(n())) %>%
mutate(sex = "Overall") %>%
ungroup() %>%
rename(age = agegrp)
metrics_by_age## # A tibble: 54 × 10
## age NHIS_9state `weighted count` `weighted count_se` `weighted percent`
## <dbl> <fct> <dbl> <dbl> <dbl>
## 1 1 Never/Never 19712179. 638174. 0.660
## 2 1 Never/Former 4913474. 318508. 0.165
## 3 1 Never/Current 2558333. 232928. 0.0857
## 4 1 Former/Never 202699. 66884. 0.00679
## 5 1 Former/Former 516129. 91968. 0.0173
## 6 1 Former/Current 914901. 134160. 0.0306
## 7 1 Current/Never 152710. 50947. 0.00511
## 8 1 Current/Former 362278. 86702. 0.0121
## 9 1 Current/Current 530331. 95138. 0.0178
## 10 2 Never/Never 24745127. 566621. 0.579
## # ℹ 44 more rows
## # ℹ 5 more variables: `weighted percent_se` <dbl>,
## # `weighted percent_low` <dbl>, `weighted percent_upp` <dbl>,
## # unweighted_count <int>, sex <chr># calculation is similar to the overall survey metrics
# but it additionally groups the data by age and sex, providing outputs for each combination of sex and age within each tobacco use status combination
metrics_by_sex_age = as_survey(svy_design) %>%
group_by(male, agegrp, NHIS_9state, .drop = FALSE) %>%
summarize("weighted count" = survey_total(),
"weighted percent" = survey_mean(vartype = c("se", "ci"), level = 0.95, proportion = TRUE),
"unweighted_count" = unweighted(n())) %>%
ungroup() %>%
rename(sex = male, age = agegrp)
metrics_by_sex_age## # A tibble: 108 × 10
## sex age NHIS_9state `weighted count` `weighted count_se`
## <dbl> <dbl> <fct> <dbl> <dbl>
## 1 0 1 Never/Never 10032119. 456161.
## 2 0 1 Never/Former 2503567. 227313.
## 3 0 1 Never/Current 1265160. 164763.
## 4 0 1 Former/Never 92241. 47705.
## 5 0 1 Former/Former 244164. 59582.
## 6 0 1 Former/Current 334887. 81472.
## 7 0 1 Current/Never 58439. 28240.
## 8 0 1 Current/Former 121394. 50683.
## 9 0 1 Current/Current 199541. 57610.
## 10 0 2 Never/Never 13863422. 420756.
## # ℹ 98 more rows
## # ℹ 5 more variables: `weighted percent` <dbl>, `weighted percent_se` <dbl>,
## # `weighted percent_low` <dbl>, `weighted percent_upp` <dbl>,
## # unweighted_count <int># combine survey metrics into a data frame
NHIS_2023_metrics<-rbind.data.frame(metrics_all,metrics_by_age,metrics_by_sex_age)
# select relevant columns and recode for readability
NHIS_2023_metrics<-NHIS_2023_metrics %>%
dplyr::select(NHIS_9state,age,sex,"weighted count","unweighted_count","weighted percent","weighted percent_se","weighted percent_low","weighted percent_upp") %>%
# filter to keep only rows where unweighted_count > 0
filter(unweighted_count>0)
NHIS_2023_metrics<-rename(NHIS_2023_metrics,cig_ecig=NHIS_9state)
NHIS_2023_metrics<-rename(NHIS_2023_metrics,"unweighted count"=unweighted_count)
# recode sex and age variables
NHIS_2023_metrics$sex[NHIS_2023_metrics$sex %in% 1]<-"Male"
NHIS_2023_metrics$sex[NHIS_2023_metrics$sex %in% 0]<-"Female"
NHIS_2023_metrics$age[NHIS_2023_metrics$age %in% 1]<-"18-24"
NHIS_2023_metrics$age[NHIS_2023_metrics$age %in% 2]<-"25-34"
NHIS_2023_metrics$age[NHIS_2023_metrics$age %in% 3]<-"35-44"
NHIS_2023_metrics$age[NHIS_2023_metrics$age %in% 4]<-"45-54"
NHIS_2023_metrics$age[NHIS_2023_metrics$age %in% 5]<-"55-64"
NHIS_2023_metrics$age[NHIS_2023_metrics$age %in% 6]<-"65+"
# export the processed data
write.csv(NHIS_2023_metrics,"NHIS_2023_metrics.csv",row.names = F)National Survey on Drug Use and Health (NSDUH)
# remove all the objects present in the workspace
rm(list=ls())
# load the data
# this sample data is a subset of NSDUH 2022
NSDUH_2022 <- read_csv("sample_data/NSDUH_2022_data.csv")
# These sample data are for instructional purposes only.
# tobacco product use status
## cigarette
## never: never smoked 100+ cigarettes in their lifetime
## current: ever smoked 100+ cigarettes in their lifetime and used at least 1 day in the past 30 days
## former: neither "current" nor "never" users (complement)
## missing: missing data
table(NSDUH_2022$UM_smkstat, useNA = "always")##
## current former missing never <NA>
## 7166 17709 59 34135 0## e-cigarette
## never: never vaped e-cigarettes
## current: ever vaped e-cigarettes and used at least 1 day in the past 30 days
## former: neither "current" nor "never" users (complement)
## missing: missing data
table(NSDUH_2022$UM_ecigstat, useNA = "always")##
## current former missing never <NA>
## 6663 10489 136 41781 0# Age: 1= <=17, 2=18-25, 3=26-34, 4=35-49, 5=50-64, 6=65+
table(NSDUH_2022$UM_agecat, useNA = "always")##
## 1 2 3 4 5 6 <NA>
## 11969 14307 9645 12588 5369 5191 0NSDUH_2022 <- NSDUH_2022[!(NSDUH_2022$UM_agecat == 1),] # exclude non-adults for this analysis
# Sex: 1=Male, 2=Female
table(NSDUH_2022$irsex, useNA = "always")##
## 1 2 <NA>
## 20852 26248 0# Filter out rows with missing value
NSDUH_2022 <- NSDUH_2022[!is.na(NSDUH_2022$UM_agecat),] # age
NSDUH_2022 <- NSDUH_2022[!is.na(NSDUH_2022$irsex),] # sex# generate a combination of tobacco product use status for cigarettes, cigars, and e-cigarettes
# creates a vector that contains all possible combinations
cig_ecig_name<-paste0(rep(c("Never","Former","Current"),each=3),"/",
rep(c("Never","Former","Current")))
# loop through all combinations and assign values:
NSDUH_2022$UM_cig_ecig_9state<-NA
m<-1
for(i in c("never","former","current")){
for(j in c("never","former","current")){
NSDUH_2022$UM_cig_ecig_9state[NSDUH_2022$UM_smkstat %in% i & NSDUH_2022$UM_ecigstat %in% j]<-cig_ecig_name[m]
m<-m+1
}
}
# For example, the loop will generate and assign the combination "Current/Never" to the UM_cig_ecig_9state column for observations where cigarettes are "Current", and e-cigarettes are "Never".
NSDUH_2022$UM_cig_ecig_9state<-factor(NSDUH_2022$UM_cig_ecig_9state,levels=unique(cig_ecig_name))
table(NSDUH_2022$UM_cig_ecig_9state,useNA = "always")##
## Never/Never Never/Former Never/Current Former/Never Former/Former
## 19765 2353 901 10106 4095
## Former/Current Current/Never Current/Former Current/Current <NA>
## 2788 2377 2507 2095 113NSDUH_2022<-NSDUH_2022[!is.na(NSDUH_2022$UM_cig_ecig_9state),]# set options for survey analysis
options(survey.lonely.psu = "adjust") # adjusts for lonely primary sampling units (PSUs) by redistributing weights
options(warn = -1) # suppress warnings
options(digits = 3) # shows 3 significant digits to display for numerical output
# Create a survey design object using the specified survey design parameters
svy_design <- svydesign(
id = ~verep, # primary sampling unit (PSU) variable
strata = ~VESTR_C, # variance stratum variable
weights = ~ANALWT2_C, # person level sampling weights
data = NSDUH_2022, # your data set
nest = TRUE # specify that the strata and PSUs are nested
)# calculate overall survey metrics by grouping the data by tobacco use status combinations
# compute weighted counts, percentages with 95% confidence intervals, and unweighted counts
# add overall indicators for sex and age
metrics_all <- as_survey(svy_design) %>% # convert the svy_design to a survey design object
# group by tobacco use status combinations, and ensures that all levels are included
group_by(UM_cig_ecig_9state, .drop = FALSE) %>%
# calculate the outputs for each group
summarize("weighted count" = survey_total(),
"weighted percent" = survey_mean(vartype = c("se","ci"), level = 0.95, proportion = TRUE),
"unweighted_count" = unweighted(n())) %>%
# indicate that the metrics are for which group, here for the overall population
mutate(sex = "Overall",age = "Overall") %>%
ungroup()
metrics_all## # A tibble: 9 × 10
## UM_cig_ecig_9state `weighted count` `weighted count_se` `weighted percent`
## <fct> <dbl> <dbl> <dbl>
## 1 Never/Never 102737571. 2425166. 0.402
## 2 Never/Former 7155393. 280089. 0.0280
## 3 Never/Current 2496607. 165234. 0.00976
## 4 Former/Never 72886624. 1533851. 0.285
## 5 Former/Former 19141320. 589288. 0.0748
## 6 Former/Current 10893035. 469528. 0.0426
## 7 Current/Never 17414043. 551115. 0.0681
## 8 Current/Former 14649985. 635456. 0.0573
## 9 Current/Current 8417375. 384883. 0.0329
## # ℹ 6 more variables: `weighted percent_se` <dbl>,
## # `weighted percent_low` <dbl>, `weighted percent_upp` <dbl>,
## # unweighted_count <int>, sex <chr>, age <chr># calculation is similar to the overall survey metrics
# but it additionally groups the data by age categories, providing outputs for each age group within each tobacco use status combination
metrics_by_age = as_survey(svy_design) %>%
group_by(UM_agecat, UM_cig_ecig_9state, .drop = FALSE) %>%
summarize("weighted count" = survey_total(),
"weighted percent" = survey_mean(vartype = c("se", "ci"), level = 0.95, proportion = TRUE),
"unweighted_count" = unweighted(n())) %>%
mutate(sex = "Overall") %>%
ungroup() %>%
rename(age = UM_agecat)
metrics_by_age## # A tibble: 45 × 10
## age UM_cig_ecig_9state `weighted count` `weighted count_se`
## <dbl> <fct> <dbl> <dbl>
## 1 2 Never/Never 16029704. 491061.
## 2 2 Never/Former 4684829. 190307.
## 3 2 Never/Current 1797715. 129438.
## 4 2 Former/Never 1124367. 73014.
## 5 2 Former/Former 3161943. 136757.
## 6 2 Former/Current 4160229. 201483.
## 7 2 Current/Never 319234. 51707.
## 8 2 Current/Former 983401. 77478.
## 9 2 Current/Current 2370537. 169481.
## 10 3 Never/Never 16068818. 611804.
## # ℹ 35 more rows
## # ℹ 6 more variables: `weighted percent` <dbl>, `weighted percent_se` <dbl>,
## # `weighted percent_low` <dbl>, `weighted percent_upp` <dbl>,
## # unweighted_count <int>, sex <chr># calculation is similar to the overall survey metrics
# but it additionally groups the data by age and sex, providing outputs for each combination of sex, age, and tobacco use status within each tobacco use status combination
metrics_by_sex_age = as_survey(svy_design) %>%
group_by(irsex, UM_agecat, UM_cig_ecig_9state, .drop = FALSE) %>%
summarize("weighted count" = survey_total(),
"weighted percent" = survey_mean(vartype = c("se", "ci"), level = 0.95, proportion = TRUE),
"unweighted_count" = unweighted(n())) %>%
ungroup() %>%
rename(sex = irsex, age = UM_agecat)
metrics_by_sex_age## # A tibble: 90 × 10
## sex age UM_cig_ecig_9state `weighted count` `weighted count_se`
## <dbl> <dbl> <fct> <dbl> <dbl>
## 1 1 2 Never/Never 7933858. 303049.
## 2 1 2 Never/Former 2079727. 126909.
## 3 1 2 Never/Current 873926. 102608.
## 4 1 2 Former/Never 606147. 58597.
## 5 1 2 Former/Former 1594545. 88770.
## 6 1 2 Former/Current 2114993. 132776.
## 7 1 2 Current/Never 218984. 45810.
## 8 1 2 Current/Former 610019. 57773.
## 9 1 2 Current/Current 1375245. 110109.
## 10 1 3 Never/Never 6813384. 384701.
## # ℹ 80 more rows
## # ℹ 5 more variables: `weighted percent` <dbl>, `weighted percent_se` <dbl>,
## # `weighted percent_low` <dbl>, `weighted percent_upp` <dbl>,
## # unweighted_count <int># combine survey metrics into a data frame
NSDUH_2022_metrics<-rbind.data.frame(metrics_all,metrics_by_age,metrics_by_sex_age)
# select relevant columns and recode for readability
NSDUH_2022_metrics<-NSDUH_2022_metrics %>%
dplyr::select(UM_cig_ecig_9state,age,sex,"weighted count","unweighted_count","weighted percent","weighted percent_se","weighted percent_low","weighted percent_upp") %>%
# filter to keep only rows where unweighted_count > 0
filter(unweighted_count>0)
NSDUH_2022_metrics<-rename(NSDUH_2022_metrics,cig_ecig=UM_cig_ecig_9state)
NSDUH_2022_metrics<-rename(NSDUH_2022_metrics,"unweighted count"=unweighted_count)
# recode sex and age variables
NSDUH_2022_metrics$sex[NSDUH_2022_metrics$sex %in% 1]<-"Male"
NSDUH_2022_metrics$sex[NSDUH_2022_metrics$sex %in% 2]<-"Female"
NSDUH_2022_metrics$age[NSDUH_2022_metrics$age %in% 2]<-"18-25"
NSDUH_2022_metrics$age[NSDUH_2022_metrics$age %in% 3]<-"26-34"
NSDUH_2022_metrics$age[NSDUH_2022_metrics$age %in% 4]<-"35-49"
NSDUH_2022_metrics$age[NSDUH_2022_metrics$age %in% 5]<-"50-64"
NSDUH_2022_metrics$age[NSDUH_2022_metrics$age %in% 6]<-"65+"
# export the processed data
write.csv(NSDUH_2022_metrics,"NSDUH_2022_metrics.csv",row.names = F)Monitoring The Future (MTF)
# remove all the objects present in the workspace
rm(list=ls())
# load the data
# this sample data is a subset of MTF 2022
MTF_2022 <- read_csv("sample_data/MTF_2022_data.csv")
# These sample data are for instructional purposes only. This instructional use has been reviewed and approved by the ICPSR Data Stewardship Policy Committee.
# tobacco product use status
## cigarette
## never: never smoked cigarettes
## current: ever smoked cigarettes and used at least 1 day in the past 30 days
## former: neither "current" nor "never" users (complement)
## missing: missing data
table(MTF_2022$cig_st, useNA = "always")##
## current former missing never <NA>
## 591 2576 1648 26623 0## e-cigarette
## never: never vaped e-cigarettes
## current: ever vaped e-cigarettes and used at least 1 day in the past 30 days
## former: neither "current" nor "never" users (complement)
## missing: missing data
table(MTF_2022$ecig_st, useNA = "always")##
## current former missing never <NA>
## 4011 4040 2186 21201 0# Sex: Male, Female
table(MTF_2022$sex, useNA = "always")##
## Female Male Unknown <NA>
## 13188 14214 4036 0# grade: 8, 10, 12
table(MTF_2022$grade, useNA = "always")##
## 8 10 12 <NA>
## 9889 11950 9599 0# Filter out rows with missing value
MTF_2022<-MTF_2022[MTF_2022$sex!="Unknown",]
table(MTF_2022$sex, useNA="always") # sex##
## Female Male <NA>
## 13188 14214 0MTF_2022<-MTF_2022[!is.na(MTF_2022$grade),]
table(MTF_2022$grade, useNA="always") #grade##
## 8 10 12 <NA>
## 8616 10311 8475 0# generate a combination of tobacco product use status for cigarettes and e-cigarettes
# creates a vector that contains all possible combinations
cig_ecig_name<-paste0(rep(c("Never","Former","Current"),each=3),"/",c("Never","Former","Current"))
# loop through all combinations and assign values:
MTF_2022$MTF_9state<-NA
k<-1
for(i in c("never","former","current")){
for(j in c("never","former","current")){
MTF_2022$MTF_9state[MTF_2022$cig_st %in% i & MTF_2022$ecig_st %in% j]<-cig_ecig_name[k]
k<-k+1
}
}
# For example, the loop will generate and assign the combination "Current/Former" to the MTF_9state column for observations where cigarettes are "Current" and e-cigarettes are "Former".
MTF_2022$MTF_9state<-factor(MTF_2022$MTF_9state, levels=unique(cig_ecig_name))
MTF_2022<-MTF_2022[!is.na(MTF_2022$MTF_9state),]
table(MTF_2022$MTF_9state,useNA = "always")##
## Never/Never Never/Former Never/Current Former/Never Former/Former
## 19376 3017 2134 407 682
## Former/Current Current/Never Current/Former Current/Current <NA>
## 1168 83 50 337 0# set options for survey analysis
options(survey.lonely.psu = "adjust") # adjusts for lonely primary sampling units (PSUs) by redistributing weights
options(warn = -1) # suppress warnings
options(digits = 3) # shows 3 significant digits to display for numerical output
# Create a survey design object using the specified survey design parameters
svy_design <- svydesign(
id = ~psu, # primary sampling unit (PSU) variable
strata = ~strata, # variance stratum variable
weights = ~weight, # person level sampling weights
data = MTF_2022, # your data set
nest = TRUE # specify that the strata and PSUs are nested
)# calculate overall survey metrics by grouping the data by tobacco use status combinations
# compute weighted counts, percentages with 95% confidence intervals, and unweighted counts
# add overall indicators for sex and age
metrics_all <- as_survey(svy_design) %>% # convert the svy_design to a survey design object
# group by tobacco use status combinations, and ensures that all levels are included
group_by(MTF_9state, .drop = FALSE) %>%
# calculate the outputs for each group
summarize("weighted count" = survey_total(),
"weighted percent" = survey_mean(vartype = c("se","ci"), level = 0.95, proportion = TRUE),
"unweighted_count" = unweighted(n())) %>%
# indicate that the metrics are for which group, here for the overall population
mutate(sex = "Overall",grade = "Overall") %>%
ungroup()
metrics_all## # A tibble: 9 × 10
## MTF_9state `weighted count` `weighted count_se` `weighted percent`
## <fct> <dbl> <dbl> <dbl>
## 1 Never/Never 19317. 3119. 0.704
## 2 Never/Former 3096. 495. 0.113
## 3 Never/Current 2178. 327. 0.0793
## 4 Former/Never 429. 90.2 0.0156
## 5 Former/Former 687. 145. 0.0250
## 6 Former/Current 1224. 215. 0.0446
## 7 Current/Never 77.2 21.0 0.00281
## 8 Current/Former 52.6 13.5 0.00191
## 9 Current/Current 394. 103. 0.0143
## # ℹ 6 more variables: `weighted percent_se` <dbl>,
## # `weighted percent_low` <dbl>, `weighted percent_upp` <dbl>,
## # unweighted_count <int>, sex <chr>, grade <chr># calculation is similar to the overall survey metrics
# but it additionally groups the data by sex categories, providing outputs for each sex group within each tobacco use status combination
metrics_by_sex <- as_survey(svy_design) %>%
group_by(sex, MTF_9state, .drop = FALSE) %>%
summarize("weighted count" = survey_total(),
"weighted percent" = survey_mean(vartype = c("se", "ci"), level = 0.95, proportion = TRUE),
"unweighted_count" = unweighted(n())) %>%
mutate(grade = "Overall") %>%
ungroup() %>%
rename(sex = sex) # rename column names if needed
metrics_by_sex## # A tibble: 18 × 10
## sex MTF_9state `weighted count` `weighted count_se` `weighted percent`
## <chr> <fct> <dbl> <dbl> <dbl>
## 1 Female Never/Never 8903. 1427. 0.674
## 2 Female Never/Former 1721. 310. 0.130
## 3 Female Never/Current 1301. 209. 0.0986
## 4 Female Former/Never 160. 32.2 0.0122
## 5 Female Former/Former 324. 61.6 0.0245
## 6 Female Former/Current 589. 107. 0.0446
## 7 Female Current/Never 18.8 6.45 0.00142
## 8 Female Current/Former 20.7 5.73 0.00157
## 9 Female Current/Curre… 164. 27.2 0.0124
## 10 Male Never/Never 10413. 1699. 0.731
## 11 Male Never/Former 1375. 190. 0.0965
## 12 Male Never/Current 876. 122. 0.0615
## 13 Male Former/Never 269. 60.2 0.0189
## 14 Male Former/Former 363. 84.8 0.0255
## 15 Male Former/Current 635. 117. 0.0446
## 16 Male Current/Never 58.4 16.9 0.00410
## 17 Male Current/Former 31.8 10.6 0.00223
## 18 Male Current/Curre… 230. 77.7 0.0161
## # ℹ 5 more variables: `weighted percent_se` <dbl>,
## # `weighted percent_low` <dbl>, `weighted percent_upp` <dbl>,
## # unweighted_count <int>, grade <chr># calculation is similar to the overall survey metrics
# but it additionally groups the data by grade and sex, providing outputs for each combination of sex and grade within each tobacco use status combination
metrics_by_sex_grade <- as_survey(svy_design) %>%
group_by(sex, grade, MTF_9state, .drop = FALSE) %>%
summarize("weighted count" = survey_total(),
"weighted percent" = survey_mean(vartype = c("se", "ci"), level = 0.95, proportion = TRUE),
"unweighted_count" = unweighted(n())) %>%
ungroup() %>%
rename(sex = sex, grade = grade) # rename column names if needed
metrics_by_sex_grade## # A tibble: 54 × 10
## sex grade MTF_9state `weighted count` `weighted count_se`
## <chr> <dbl> <fct> <dbl> <dbl>
## 1 Female 8 Never/Never 3299. 572.
## 2 Female 8 Never/Former 384. 118.
## 3 Female 8 Never/Current 222. 80.1
## 4 Female 8 Former/Never 43.4 10.1
## 5 Female 8 Former/Former 64.2 20.6
## 6 Female 8 Former/Current 97.4 31.4
## 7 Female 8 Current/Never 0.716 0.716
## 8 Female 8 Current/Former 5.13 2.31
## 9 Female 8 Current/Current 18.0 8.37
## 10 Female 10 Never/Never 3226. 559.
## # ℹ 44 more rows
## # ℹ 5 more variables: `weighted percent` <dbl>, `weighted percent_se` <dbl>,
## # `weighted percent_low` <dbl>, `weighted percent_upp` <dbl>,
## # unweighted_count <int># combine survey metrics into a data frame
MTF_2022_metrics<-rbind.data.frame(metrics_all,metrics_by_sex,metrics_by_sex_grade)
# select relevant columns and recode for readability
MTF_2022_metrics<-MTF_2022_metrics %>%
dplyr::select(MTF_9state,grade,sex,"weighted count","unweighted_count","weighted percent","weighted percent_se","weighted percent_low","weighted percent_upp") %>%
# filter to keep only rows where unweighted_count > 0
filter(unweighted_count>0)
MTF_2022_metrics<-rename(MTF_2022_metrics,cig_ecig=MTF_9state)
MTF_2022_metrics<-rename(MTF_2022_metrics,"unweighted count"=unweighted_count)
# recode sex and grade variables if needed
# export the processed data
write.csv(MTF_2022_metrics,"MTF_2022_metrics.csv",row.names = F)National Youth Tobacco Survey (NYTS)
# remove all the objects present in the workspace
rm(list=ls())
# load the data
# this sample data is a subset of NYTS 2023
NYTS_2023 <- read_csv("sample_data/NYTS_2023_data.csv")
# These sample data are for instructional purposes only.
# tobacco product use status
## cigarette
## never: never smoked 100+ cigarettes in their lifetime
## current: ever smoked 100+ cigarettes in their lifetime and used at least 1 day in the past 30 days
## former: neither "current" nor "never" users (complement)
## missing: missing data
table(NYTS_2023$cigt, useNA = "always")##
## Current Former Never <NA>
## 80 42 20006 0## e-cigarette
## never: never vaped e-cigarettes
## current: ever vaped e-cigarettes and used at least 1 day in the past 30 days
## former: neither "current" nor "never" users (complement)
## missing: missing data
table(NYTS_2023$ecig, useNA = "always")##
## Current Former Never <NA>
## 1382 1631 17115 0# Age: 12-14, 15-17, 18+
table(NYTS_2023$age_cat, useNA = "always")##
## 12-14 15-17 18+ <NA>
## 10517 8506 1105 0# Sex: Male, Female
table(NYTS_2023$um_sex_lab, useNA = "always")##
## Female Male <NA>
## 9859 10269 0# Filter out rows with missing value
NYTS_2023 <- NYTS_2023[!is.na(NYTS_2023$age_cat),] # age
NYTS_2023 <- NYTS_2023[!is.na(NYTS_2023$um_sex_lab),] # sex# generate a combination of tobacco product use status for cigarettes and e-cigarettes
# creates a vector that contains all possible combinations
cig_ecig_name<-paste0(rep(c("Never","Former","Current"),each=3),"/",c("Never","Former","Current"))
# loop through all combinations and assign values:
NYTS_2023$NYTS_9state<-NA
k<-1
for(i in c("Never","Former","Current")){
for(j in c("Never","Former","Current")){
NYTS_2023$NYTS_9state[NYTS_2023$cigt %in% i & NYTS_2023$ecig %in% j]<-cig_ecig_name[k]
k<-k+1
}
}
# For example, the loop will generate and assign the combination "Current/Former" to the NYTS_9state column for observations where cigarettes are "Current" and e-cigarettes are "Former".
NYTS_2023$NYTS_9state<-factor(NYTS_2023$NYTS_9state, levels=unique(cig_ecig_name))
NYTS_2023<-NYTS_2023[!is.na(NYTS_2023$NYTS_9state),]
table(NYTS_2023$NYTS_9state,useNA = "always")##
## Never/Never Never/Former Never/Current Former/Never Former/Former
## 17092 1620 1294 6 9
## Former/Current Current/Never Current/Former Current/Current <NA>
## 27 17 2 61 0# set options for survey analysis
options(survey.lonely.psu = "adjust") # adjusts for lonely primary sampling units (PSUs) by redistributing weights
options(warn = -1) # suppress warnings
options(digits = 3) # shows 3 significant digits to display for numerical output
# Create a survey design object using the specified survey design parameters
svy_design <- svydesign(
id = ~psu, # primary sampling unit (PSU) variable
strata = ~strata, # variance stratum variable
weights = ~weights, # person-level sampling weights
data = NYTS_2023, # your data set
nest = TRUE # specify that the strata and PSUs are nested
)# calculate overall survey metrics by grouping the data by tobacco use status combinations
# compute weighted counts, percentages with 95% confidence intervals, and unweighted counts
# add overall indicators for sex and age
metrics_all <- as_survey(svy_design) %>% # convert the svy_design to a survey design object
# group by tobacco use status combinations, and ensures that all levels are included
group_by(NYTS_9state, .drop = FALSE) %>%
# calculate the outputs for each group
summarize("weighted count" = survey_total(),
"weighted percent" = survey_mean(vartype = c("se","ci"), level = 0.95, proportion = TRUE),
"unweighted_count" = unweighted(n())) %>%
# indicate that the metrics are for which group, here for the overall population
mutate(sex = "Overall",age = "Overall") %>%
ungroup()
metrics_all## # A tibble: 9 × 10
## NYTS_9state `weighted count` `weighted count_se` `weighted percent`
## <fct> <dbl> <dbl> <dbl>
## 1 Never/Never 21811551. 1470758. 0.836
## 2 Never/Former 2320582. 226331. 0.0889
## 3 Never/Current 1809811. 160046. 0.0694
## 4 Former/Never 3467. 1861. 0.000133
## 5 Former/Former 10922. 4598. 0.000419
## 6 Former/Current 35906. 10750. 0.00138
## 7 Current/Never 32199. 15729. 0.00123
## 8 Current/Former 3330. 2587. 0.000128
## 9 Current/Current 60935. 13951. 0.00234
## # ℹ 6 more variables: `weighted percent_se` <dbl>,
## # `weighted percent_low` <dbl>, `weighted percent_upp` <dbl>,
## # unweighted_count <int>, sex <chr>, age <chr># calculation is similar to the overall survey metrics
# but it additionally groups the data by age categories, providing outputs for each age group within each tobacco use status combination
metrics_by_age <- as_survey(svy_design) %>%
group_by(age_cat, NYTS_9state, .drop = FALSE) %>%
summarize("weighted count" = survey_total(),
"weighted percent" = survey_mean(vartype = c("se", "ci"), level = 0.95, proportion = TRUE),
"unweighted_count" = unweighted(n())) %>%
mutate(sex = "Overall") %>%
ungroup() %>%
rename(age = age_cat)
metrics_by_age## # A tibble: 27 × 10
## age NYTS_9state `weighted count` `weighted count_se` `weighted percent`
## <chr> <fct> <dbl> <dbl> <dbl>
## 1 12-14 Never/Never 10630294. 1017354. 0.903
## 2 12-14 Never/Former 571979. 62000. 0.0486
## 3 12-14 Never/Current 523608. 77380. 0.0445
## 4 12-14 Former/Never 1774. 1250. 0.000151
## 5 12-14 Former/Former 1710. 1193. 0.000145
## 6 12-14 Former/Current 9183. 7727. 0.000780
## 7 12-14 Current/Never 16172. 13962. 0.00137
## 8 12-14 Current/Former 3330. 2587. 0.000283
## 9 12-14 Current/Current 16622. 6684. 0.00141
## 10 15-17 Never/Never 9439850. 995333. 0.783
## # ℹ 17 more rows
## # ℹ 5 more variables: `weighted percent_se` <dbl>,
## # `weighted percent_low` <dbl>, `weighted percent_upp` <dbl>,
## # unweighted_count <int>, sex <chr># calculation is similar to the overall survey metrics
# but it additionally groups the data by age and sex, providing outputs for each combination of sex and age within each tobacco use status combination
metrics_by_sex_age <- as_survey(svy_design) %>%
group_by(um_sex_lab, age_cat, NYTS_9state, .drop = FALSE) %>%
summarize("weighted count" = survey_total(),
"weighted percent" = survey_mean(vartype = c("se", "ci"), level = 0.95, proportion = TRUE),
"unweighted_count" = unweighted(n())) %>%
ungroup() %>%
rename(sex = um_sex_lab, age = age_cat)
metrics_by_sex_age## # A tibble: 54 × 10
## sex age NYTS_9state `weighted count` `weighted count_se`
## <chr> <chr> <fct> <dbl> <dbl>
## 1 Female 12-14 Never/Never 5212285. 537237.
## 2 Female 12-14 Never/Former 303541. 42140.
## 3 Female 12-14 Never/Current 300967. 46213.
## 4 Female 12-14 Former/Never 0 0
## 5 Female 12-14 Former/Former 802. 775.
## 6 Female 12-14 Former/Current 1630. 1630.
## 7 Female 12-14 Current/Never 1646. 1122.
## 8 Female 12-14 Current/Former 0 0
## 9 Female 12-14 Current/Current 6258. 4873.
## 10 Female 15-17 Never/Never 4512586. 537851.
## # ℹ 44 more rows
## # ℹ 5 more variables: `weighted percent` <dbl>, `weighted percent_se` <dbl>,
## # `weighted percent_low` <dbl>, `weighted percent_upp` <dbl>,
## # unweighted_count <int># combine survey metrics into a data frame
NYTS_2023_metrics<-rbind.data.frame(metrics_all,metrics_by_age,metrics_by_sex_age)
# select relevant columns and recode for readability
NYTS_2023_metrics<-NYTS_2023_metrics %>%
dplyr::select(NYTS_9state,age,sex,"weighted count","unweighted_count","weighted percent","weighted percent_se","weighted percent_low","weighted percent_upp") %>%
# filter to keep only rows where unweighted_count > 0
filter(unweighted_count>0)
NYTS_2023_metrics<-rename(NYTS_2023_metrics,cig_ecig=NYTS_9state)
NYTS_2023_metrics<-rename(NYTS_2023_metrics,"unweighted count"=unweighted_count)
# recode sex and age variables if needed
# export the processed data
write.csv(NYTS_2023_metrics,"NYTS_2023_metrics.csv",row.names = F)The Tobacco Use Supplement to the Current Population Survey (TUS-CPS)
rm(list=ls()) # remove all the objects present in the workspace
# load the data
# this data set is a subset of TUSCPS 2022
TUSCPS_2022 <- read_csv("sample_data/TUSCPS_2022rep_data.csv")
# tobacco product use status
## cigarette
## 1: current: ever smoked 100+ cigarettes in their lifetime and used at least 1 day in the past 30 days
## 2: former: neither "current" nor "never" users (complement)
## 3: never: never smoked 100+ cigarettes in their lifetime
table(TUSCPS_2022$cig_status, useNA = "always")##
## current former missing never <NA>
## 87 83 3 249 0## e-cigarette
## 1: current: ever vaped e-cigarettes and used at least 1 day in the past 30 days
## 2: former: neither "current" nor "never" users (complement)
## 3: never: never vaped e-cigarettes
table(TUSCPS_2022$ecig_status, useNA = "always")##
## current former missing never <NA>
## 14 50 10 348 0# Age: 1=18-24, 2=25-34, 3=35-44, 4=45-54, 5=55-64, 6=65+
table(TUSCPS_2022$Age, useNA = "always")##
## 1 2 3 4 5 6 <NA>
## 25 90 79 71 60 97 0# Sex:1=Male, 2=Female
table(TUSCPS_2022$Sex, useNA = "always")##
## 1 2 <NA>
## 166 256 0# Filter out rows with missing value
TUSCPS_2022 <- TUSCPS_2022[!is.na(TUSCPS_2022$Age),] # age
TUSCPS_2022 <- TUSCPS_2022[!is.na(TUSCPS_2022$Sex),] # sex# generate a combination of tobacco product use status for cigarettes and e-cigarettes
# creates a vector that contains all possible combinations
cig_ecig_name<-paste0(rep(c("Never","Former","Current"),each=3),"/",c("Never","Former","Current"))
# loop through all combinations and assign values:
TUSCPS_2022$TUSCPS_9state<-NA
k<-1
for(i in c("never","former","current")){
for(j in c("never","former","current")){
TUSCPS_2022$TUSCPS_9state[TUSCPS_2022$cig_status %in% i & TUSCPS_2022$ecig_status %in% j]<-cig_ecig_name[k]
k<-k+1
}
}
# For example, the loop will generate and assign the combination "Current/Former" to the TUSCPS_9state column for observations where cigarettes are "Current" and e-cigarettes are "Former"
TUSCPS_2022$TUSCPS_9state<-factor(TUSCPS_2022$TUSCPS_9state, levels=unique(cig_ecig_name))
TUSCPS_2022 <- TUSCPS_2022[!is.na(TUSCPS_2022$TUSCPS_9state),]
table(TUSCPS_2022$TUSCPS_9state,useNA = "always")##
## Never/Never Never/Former Never/Current Former/Never Former/Former
## 231 13 1 65 13
## Former/Current Current/Never Current/Former Current/Current <NA>
## 3 51 24 10 0# set options for survey analysis
options(survey.replicates.mse = TRUE) # compute the mean square error for replicate weights
options(warn = -1) # suppress warnings
options(digits = 3) # shows 3 significant digits to display for numerical output
# create a survey design object using replicate weights
repweight <- TUSCPS_2022[, paste0("pwsrwgt", 1:160)]
svy_design <- svrepdesign(
weights = ~PWSRWGT, # sampling weights
repweights = repweight, # replicate weights for TUS-CPS
type = "BRR", # replication method
fay.rho = 0.5, # set the rho parameter for Fay's method
data = TUSCPS_2022 # your data set
)# calculate overall survey metrics by grouping the data by tobacco use status combinations
# compute weighted counts, percentages with 95% confidence intervals, and unweighted counts
# add overall indicators for sex and age
metrics_all <- as_survey_rep(svy_design) %>% # use as_survey_rep for replicate weights
# group by tobacco use status combinations, and ensures that all levels are included
group_by(TUSCPS_9state, .drop = FALSE) %>%
# calculate the outputs for each group
summarize("weighted count" = survey_total(),
"weighted percent" = survey_mean(vartype = c("se","ci"), level = 0.95, proportion = TRUE),
"unweighted_count" = unweighted(n())) %>%
# indicate that the metrics are for which group, here for the overall population
mutate(sex = "Overall",age = "Overall") %>%
ungroup()
metrics_all## # A tibble: 9 × 10
## TUSCPS_9state `weighted count` `weighted count_se` `weighted percent`
## <fct> <dbl> <dbl> <dbl>
## 1 Never/Never 1486492. 14859464678. 0.574
## 2 Never/Former 88446. 878277928. 0.0341
## 3 Never/Current 14738. 163031941. 0.00569
## 4 Former/Never 376635. 3771480956. 0.145
## 5 Former/Former 70903. 714887534. 0.0274
## 6 Former/Current 38429. 407044512. 0.0148
## 7 Current/Never 257866. 2591164913. 0.0995
## 8 Current/Former 165959. 1663405998. 0.0640
## 9 Current/Current 92278. 930598832. 0.0356
## # ℹ 6 more variables: `weighted percent_se` <dbl>,
## # `weighted percent_low` <dbl>, `weighted percent_upp` <dbl>,
## # unweighted_count <int>, sex <chr>, age <chr># calculation is similar to the overall survey metrics
# but it additionally groups the data by age categories, providing outputs for each age group within each tobacco use status combination
metrics_by_age <- as_survey_rep(svy_design) %>%
group_by(Age, TUSCPS_9state, .drop = FALSE) %>%
summarize("weighted count" = survey_total(),
"weighted percent" = survey_mean(vartype = c("se", "ci"), level = 0.95, proportion = TRUE),
"unweighted_count" = unweighted(n())) %>%
mutate(sex = "Overall") %>%
ungroup() %>%
rename(age = Age)
metrics_by_age## # A tibble: 54 × 10
## age TUSCPS_9state `weighted count` `weighted count_se` `weighted percent`
## <dbl> <fct> <dbl> <dbl> <dbl>
## 1 1 Never/Never 117605. 1157511717. 5.07e- 1
## 2 1 Never/Former 20664. 221237207. 8.91e- 2
## 3 1 Never/Current 14738. 163031941. 6.35e- 2
## 4 1 Former/Never 7872. 85489628. 3.39e- 2
## 5 1 Former/Former 1897. 21471531. 8.18e- 3
## 6 1 Former/Current 0 0 5.89e-12
## 7 1 Current/Never 29331. 304261991. 1.26e- 1
## 8 1 Current/Former 9962. 107204862. 4.30e- 2
## 9 1 Current/Current 29856. 320274441. 1.29e- 1
## 10 2 Never/Never 408645. 4093026080. 6.36e- 1
## # ℹ 44 more rows
## # ℹ 5 more variables: `weighted percent_se` <dbl>,
## # `weighted percent_low` <dbl>, `weighted percent_upp` <dbl>,
## # unweighted_count <int>, sex <chr># calculation is similar to the overall survey metrics
# but it additionally groups the data by age and sex, providing outputs for each combination of sex and age within each tobacco use status combination
metrics_by_sex_age <- as_survey_rep(svy_design) %>%
group_by(Sex, Age, TUSCPS_9state, .drop = FALSE) %>%
summarize("weighted count" = survey_total(),
"weighted percent" = survey_mean(vartype = c("se", "ci"), level = 0.95, proportion = TRUE),
"unweighted_count" = unweighted(n())) %>%
ungroup() %>%
rename(sex = Sex, age = Age)
metrics_by_sex_age## # A tibble: 108 × 10
## sex age TUSCPS_9state `weighted count` `weighted count_se`
## <dbl> <dbl> <fct> <dbl> <dbl>
## 1 1 1 Never/Never 85087. 839252518.
## 2 1 1 Never/Former 0 0
## 3 1 1 Never/Current 0 0
## 4 1 1 Former/Never 4557. 51535184.
## 5 1 1 Former/Former 0 0
## 6 1 1 Former/Current 0 0
## 7 1 1 Current/Never 5580. 60766735.
## 8 1 1 Current/Former 9962. 107204862.
## 9 1 1 Current/Current 28451. 307734784.
## 10 1 2 Never/Never 162327. 1624978202.
## # ℹ 98 more rows
## # ℹ 5 more variables: `weighted percent` <dbl>, `weighted percent_se` <dbl>,
## # `weighted percent_low` <dbl>, `weighted percent_upp` <dbl>,
## # unweighted_count <int># combine survey metrics into a data frame
data <- rbind.data.frame(metrics_all, metrics_by_age, metrics_by_sex_age)
# select relevant columns and recode for readability
TUSCPS_2022_metrics <- data %>%
dplyr::select(TUSCPS_9state, age, sex,
"weighted count", "unweighted_count", "weighted percent",
"weighted percent_se", "weighted percent_low", "weighted percent_upp") %>%
# filters out rows with zero unweighted counts
filter(unweighted_count > 0)
# rename columns
TUSCPS_2022_metrics <- rename(TUSCPS_2022_metrics, cig_ecig = TUSCPS_9state)
TUSCPS_2022_metrics <- rename(TUSCPS_2022_metrics, "unweighted count" = unweighted_count)
# recode sex and age variables
TUSCPS_2022_metrics$sex[TUSCPS_2022_metrics$sex %in% 1] <- "Male"
TUSCPS_2022_metrics$sex[TUSCPS_2022_metrics$sex %in% 2] <- "Female"
TUSCPS_2022_metrics$age[TUSCPS_2022_metrics$age %in% 1] <- "18-24"
TUSCPS_2022_metrics$age[TUSCPS_2022_metrics$age %in% 2] <- "25-34"
TUSCPS_2022_metrics$age[TUSCPS_2022_metrics$age %in% 3] <- "35-44"
TUSCPS_2022_metrics$age[TUSCPS_2022_metrics$age %in% 4] <- "45-54"
TUSCPS_2022_metrics$age[TUSCPS_2022_metrics$age %in% 5] <- "55-64"
TUSCPS_2022_metrics$age[TUSCPS_2022_metrics$age %in% 6] <- "65+"
# export the processed data
write.csv(TUSCPS_2022_metrics,"TUSCPS_2022_metrics.csv",row.names = F)Youth Risk Behavior Surveillance System
(YRBSS)
Youth Risk Behavior Survey (YRBS)
# remove all the objects present in the workspace
rm(list=ls())
# load the data
# this sample data is a subset of YRBS 2021
YRBS_2021 <- read_csv("sample_data/YRBS_2021_data.csv")
# These sample data are for instructional purposes only.
# tobacco product use status
## cigarette
## never: never smoked cigarettes
## current: ever smoked cigarettes and used at least 1 day in the past 30 days
## former: neither "current" nor "never" users (complement)
## missing: missing data
table(YRBS_2021$cig_st, useNA = "always")##
## current former missing never <NA>
## 646 1735 3536 11156 0## e-cigarette
## never: never vaped e-cigarettes
## current: ever vaped e-cigarettes and used at least 1 day in the past 30 days
## former: neither "current" nor "never" users (complement)
## missing: missing data
table(YRBS_2021$ecig_st, useNA = "always")##
## current former missing never <NA>
## 2892 2262 1243 10676 0# Age: ~12, 13, 14, 15, 16, 17, 18+, Missing
table(YRBS_2021$Q1, useNA = "always")##
## 12 years old or younger 13 years old 14 years old
## 29 61 3391
## 15 years old 16 years old 17 years old
## 4410 4257 3894
## 18 years old or older Missing <NA>
## 1011 20 0# Sex: Male, Female, Missing
table(YRBS_2021$sex, useNA = "always")##
## Female Male Missing <NA>
## 8128 8765 180 0# Filter out rows with missing value
YRBS_2021 <- YRBS_2021[YRBS_2021$Q1 != "Missing",] # age
YRBS_2021 <- YRBS_2021[!is.na(YRBS_2021$Q1),]
table(YRBS_2021$Q1, useNA="always")##
## 12 years old or younger 13 years old 14 years old
## 29 61 3391
## 15 years old 16 years old 17 years old
## 4410 4257 3894
## 18 years old or older <NA>
## 1011 0YRBS_2021 <- YRBS_2021[YRBS_2021$sex != "Missing",] # sex
YRBS_2021 <- YRBS_2021[!is.na(YRBS_2021$sex),]
table(YRBS_2021$sex, useNA="always")##
## Female Male <NA>
## 8121 8755 0# generate a combination of tobacco product use status for cigarettes and e-cigarettes
# creates a vector that contains all possible combinations
cig_ecig_name<-paste0(rep(c("Never","Former","Current"),each=3),"/",c("Never","Former","Current"))
# loop through all combinations and assign values:
YRBS_2021$YRBS_9state<-NA
k<-1
for(i in c("never","former","current")){
for(j in c("never","former","current")){
YRBS_2021$YRBS_9state[YRBS_2021$cig_st %in% i & YRBS_2021$ecig_st %in% j]<-cig_ecig_name[k]
k<-k+1
}
}
# For example, the loop will generate and assign the combination "Current/Former" to the YRBS_9state column for observations where cigarettes are "Current" and e-cigarettes are "Former".
YRBS_2021$YRBS_9state<-factor(YRBS_2021$YRBS_9state, levels=unique(cig_ecig_name))
YRBS_2021<-YRBS_2021[!is.na(YRBS_2021$YRBS_9state),]
table(YRBS_2021$YRBS_9state,useNA = "always")##
## Never/Never Never/Former Never/Current Former/Never Former/Former
## 8329 1203 956 246 511
## Former/Current Current/Never Current/Former Current/Current <NA>
## 754 28 25 538 0# set options for survey analysis
options(survey.lonely.psu = "adjust") # adjusts for lonely primary sampling units (PSUs) by redistributing weights
options(warn = -1) # suppress warnings
options(digits = 3) # shows 3 significant digits to display for numerical output
# Create a survey design object using the specified survey design parameters
svy_design <- svydesign(
id = ~PSU, # primary sampling unit (PSU) variable
strata = ~STRATUM, # variance stratum variable
weights = ~WEIGHT, # sampling weights
data = YRBS_2021, # your data set
nest = TRUE # specify that the strata and PSUs are nested
)# calculate overall survey metrics by grouping the data by tobacco use status combinations
# compute weighted counts, percentages with 95% confidence intervals, and unweighted counts
# add overall indicators for sex and age
metrics_all <- as_survey(svy_design) %>% # convert the svy_design to a survey design object
# group by tobacco use status combinations, and ensures that all levels are included
group_by(YRBS_9state, .drop = FALSE) %>%
# calculate the outputs for each group
summarize("weighted count" = survey_total(),
"weighted percent" = survey_mean(vartype = c("se","ci"), level = 0.95, proportion = TRUE),
"unweighted_count" = unweighted(n())) %>%
# indicate that the metrics are for which group, here for the overall population
mutate(sex = "Overall",age = "Overall") %>%
ungroup()
metrics_all## # A tibble: 9 × 10
## YRBS_9state `weighted count` `weighted count_se` `weighted percent`
## <fct> <dbl> <dbl> <dbl>
## 1 Never/Never 8931. 745. 0.658
## 2 Never/Former 1335. 128. 0.0983
## 3 Never/Current 1118. 118. 0.0823
## 4 Former/Never 234. 32.8 0.0172
## 5 Former/Former 546. 42.6 0.0402
## 6 Former/Current 828. 79.0 0.0610
## 7 Current/Never 27.0 8.50 0.00199
## 8 Current/Former 34.3 7.10 0.00252
## 9 Current/Current 528. 58.6 0.0389
## # ℹ 6 more variables: `weighted percent_se` <dbl>,
## # `weighted percent_low` <dbl>, `weighted percent_upp` <dbl>,
## # unweighted_count <int>, sex <chr>, age <chr># calculation is similar to the overall survey metrics
# but it additionally groups the data by age categories, providing outputs for each age group within each tobacco use status combination
metrics_by_age <- as_survey(svy_design) %>%
group_by(Q1, YRBS_9state, .drop = FALSE) %>%
summarize("weighted count" = survey_total(),
"weighted percent" = survey_mean(vartype = c("se", "ci"), level = 0.95, proportion = TRUE),
"unweighted_count" = unweighted(n())) %>%
mutate(sex = "Overall") %>%
ungroup() %>%
rename(age = Q1)
metrics_by_age## # A tibble: 63 × 10
## age YRBS_9state `weighted count` `weighted count_se` `weighted percent`
## <chr> <fct> <dbl> <dbl> <dbl>
## 1 12 years… Never/Never 2.87 1.68 5.53e- 1
## 2 12 years… Never/Form… 0 0 2.40e-11
## 3 12 years… Never/Curr… 0 0 2.40e-11
## 4 12 years… Former/Nev… 0 0 2.40e-11
## 5 12 years… Former/For… 0.410 0.410 7.90e- 2
## 6 12 years… Former/Cur… 0 0 2.40e-11
## 7 12 years… Current/Ne… 0 0 2.40e-11
## 8 12 years… Current/Fo… 0 0 2.40e-11
## 9 12 years… Current/Cu… 1.91 1.13 3.68e- 1
## 10 13 years… Never/Never 36.5 8.16 8.15e- 1
## # ℹ 53 more rows
## # ℹ 5 more variables: `weighted percent_se` <dbl>,
## # `weighted percent_low` <dbl>, `weighted percent_upp` <dbl>,
## # unweighted_count <int>, sex <chr># calculation is similar to the overall survey metrics
# but it additionally groups the data by age and sex, providing outputs for each combination of sex and age within each tobacco use status combination
metrics_by_sex_age <- as_survey(svy_design) %>%
group_by(sex, Q1, YRBS_9state, .drop = FALSE) %>%
summarize("weighted count" = survey_total(),
"weighted percent" = survey_mean(vartype = c("se", "ci"), level = 0.95, proportion = TRUE),
"unweighted_count" = unweighted(n())) %>%
ungroup() %>%
rename(sex = sex, age = Q1)
metrics_by_sex_age## # A tibble: 126 × 10
## sex age YRBS_9state `weighted count` `weighted count_se`
## <chr> <chr> <fct> <dbl> <dbl>
## 1 Female 12 years old or youn… Never/Never 0 0
## 2 Female 12 years old or youn… Never/Form… 0 0
## 3 Female 12 years old or youn… Never/Curr… 0 0
## 4 Female 12 years old or youn… Former/Nev… 0 0
## 5 Female 12 years old or youn… Former/For… 0.410 0.410
## 6 Female 12 years old or youn… Former/Cur… 0 0
## 7 Female 12 years old or youn… Current/Ne… 0 0
## 8 Female 12 years old or youn… Current/Fo… 0 0
## 9 Female 12 years old or youn… Current/Cu… 1.16 1.05
## 10 Female 13 years old Never/Never 23.5 6.68
## # ℹ 116 more rows
## # ℹ 5 more variables: `weighted percent` <dbl>, `weighted percent_se` <dbl>,
## # `weighted percent_low` <dbl>, `weighted percent_upp` <dbl>,
## # unweighted_count <int># combine survey metrics into a data frame
YRBS_2021_metrics<-rbind.data.frame(metrics_all,metrics_by_age,metrics_by_sex_age)
# select relevant columns and recode for readability
YRBS_2021_metrics<-YRBS_2021_metrics %>%
dplyr::select(YRBS_9state,age,sex,"weighted count","unweighted_count","weighted percent","weighted percent_se","weighted percent_low","weighted percent_upp") %>%
# filter to keep only rows where unweighted_count > 0
filter(unweighted_count>0)
YRBS_2021_metrics<-rename(YRBS_2021_metrics,cig_ecig=YRBS_9state)
YRBS_2021_metrics<-rename(YRBS_2021_metrics,"unweighted count"=unweighted_count)
# recode sex and age variables if needed
# export the processed data
write.csv(YRBS_2021_metrics,"YRBS_2021_metrics.csv",row.names = F)