# Sample Code

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:

1. Prepare Data: Clean and filter the data and ensure it is ready for analysis.
2. Calculate Survey Metrics: Set up a survey design object with weights and configure options for analysis. Compute weighted and unweighted metrics for tobacco use status combinations, or by other category variables.
3. Export Data: Clean and format the data, and then export the processed data to a CSV file.

Each step is detailed with specific code snippets to achieve the desired analysis.

# 1 PATH

Population Assessment of Tobacco and Health (PATH) Study

## 1.1 Prepare data

    # 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             240

    PATH_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

## 1.2 Calculate survey metrics

    # 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
      )

### 1.2.1 Overall survey metrics

    # 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>

### 1.2.2 Grouped by age

    # 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>

### 1.2.3 Grouped by age and sex

    # 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>

## 1.3 Process Data

    # 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)

# 2 NHIS

National Health Interview Survey (NHIS)

## 2.1 Prepare data

    # 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

## 2.2 Calculate survey metrics

    # 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
    )

### 2.2.1 Overall survey metrics

    # 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>

### 2.2.2 Grouped by age

    # 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>

### 2.2.3 Grouped by age and sex

    # 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>

## 2.3 Process Data

    # 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)

# 3 NSDUH

National Survey on Drug Use and Health (NSDUH)

## 3.1 Prepare data

    # 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     0

    NSDUH_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             113

    NSDUH_2022<-NSDUH_2022[!is.na(NSDUH_2022$UM_cig_ecig_9state),]

## 3.2 Calculate survey metrics

    # 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
    )

### 3.2.1 Overall survey metrics

    # 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>

### 3.2.2 Grouped by age

    # 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>

### 3.2.3 Grouped by age and sex

    # 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>

## 3.3 Process Data

    # 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)

# 4 MTF

Monitoring The Future (MTF)

## 4.1 Prepare data

    # 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      0

    MTF_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

## 4.2 Calculate survey metrics

    # 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
    )

### 4.2.1 Overall survey metrics

    # 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>

### 4.2.2 Grouped by sex

    # 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>

### 4.2.3 Grouped by grade and sex

    # 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>

## 4.3 Process Data

    # 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)

# 5 NYTS

National Youth Tobacco Survey (NYTS)

## 5.1 Prepare data

    # 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

## 5.2 Calculate survey metrics

    # 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
    )

### 5.2.1 Overall survey metrics

    # 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>

### 5.2.2 Grouped by age

    # 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>

### 5.2.3 Grouped by age and sex

    # 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>

## 5.3 Process Data

    # 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)

# 6 TUS-CPS

The Tobacco Use Supplement to the Current Population Survey (TUS-CPS)

## 6.1 Prepare data

    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

## 6.2 Calculate survey metrics

    # 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
      )

### 6.2.1 Overall survey metrics

    # 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>

### 6.2.2 Grouped by age

    # 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>

### 6.2.3 Grouped by age and sex

    # 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>

## 6.3 Process Data

    # 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)

# 7 YRBSS

Youth Risk Behavior Surveillance System (YRBSS)
Youth Risk Behavior Survey (YRBS)

## 7.1 Prepare data

    # 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                       0

    YRBS_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

## 7.2 Calculate survey metrics

    # 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
    )

### 7.2.1 Overall survey metrics

    # 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>

### 7.2.2 Grouped by age

    # 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>

### 7.2.3 Grouped by age and sex

    # 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>

## 7.3 Process Data

    # 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)