---
title: "Sample Code"
Author: "Center for the Assessment of Tobacco Regulations, Data Analysis and Modeling Core"
Version: "2026_07_06"
output: 
  html_document:
    toc: true
    toc_float: true
    number_sections: true
    theme: united
    highlight: textmate
    fig_caption: yes
---

```{=html}
<style>
p.caption {
  font-size: 2em;
}
</style>
```

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = F, fig.align = "center", message = F, warning = F)
rm(list=ls())    # remove all the objects present in the workspace
  
library("knitr")
library(survey)
library(srvyr)
library(tidyverse)
```

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.

</br>

# PATH
**Population Assessment of Tobacco and Health (PATH) Study**

## Prepare data

```{r, echo = T}
# 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")

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

# age: 1: 18-24, 2: 25-34, 3: 35-54, 4: 55+
table(PATH_wave1$UM_W1_A_agecat, useNA = "always")

# sex: 0: Female, 1: Male
table(PATH_wave1$UM_W1_A_male_imp, useNA = "always")

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

```{r, echo = T}
# 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")

PATH_wave1<-PATH_wave1[!is.na(PATH_wave1$UM_W1_cig_ecig_9state),]
table(PATH_wave1$UM_W1_cig_ecig_9state,useNA = "always")
```

## Calculate survey metrics

```{r, echo = T}
# 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
  )
```

### Overall survey metrics

```{r, echo = T}
# 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
```

### Grouped by age

```{r, echo = T}
# 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
```

### Grouped by age and sex

```{r,echo=T}
# 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
```

## Process Data

```{r,echo=T}
# 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)
```

# NHIS
**National Health Interview Survey (NHIS)**

## Prepare data

```{r, echo = T}
# 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")

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

# Age: 1=18-24, 2=25-34, 3=35-44, 4=45-54, 5=55-64, 6=65+
table(NHIS_2023$agegrp, useNA = "always")

# Sex: 0=Female, 1=Male
table(NHIS_2023$male, useNA = "always")

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

```{r, echo = T}
# 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")

```

## Calculate survey metrics

```{r, echo = T}
# 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
)
```

### Overall survey metrics

```{r, echo = T}
# 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
```

### Grouped by age

```{r, echo = T}
# 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
```

### Grouped by age and sex

```{r,echo=T}
# 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
```

## Process Data

```{r,echo=T}
# 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)
```

# NSDUH
**National Survey on Drug Use and Health (NSDUH)**

## Prepare data

```{r, echo = T}
# 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")

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

# Age: 1= <=17, 2=18-25, 3=26-34, 4=35-49, 5=50-64, 6=65+
table(NSDUH_2022$UM_agecat, useNA = "always")
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")

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

```{r, echo = T}
# 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")

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

## Calculate survey metrics

```{r, echo = T}
# 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
)
```

### Overall survey metrics

```{r, echo = T}
# 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
```

### Grouped by age

```{r, echo = T}
# 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
```

### Grouped by age and sex

```{r,echo=T}
# 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
```

## Process Data

```{r,echo=T}
# 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)
```

# MTF
**Monitoring The Future (MTF)**

## Prepare data

```{r, echo = T}
# 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")

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

# Sex: Male, Female
table(MTF_2022$sex, useNA = "always")

# grade: 8, 10, 12
table(MTF_2022$grade, useNA = "always")

# Filter out rows with missing value
MTF_2022<-MTF_2022[MTF_2022$sex!="Unknown",]
table(MTF_2022$sex, useNA="always")     # sex

MTF_2022<-MTF_2022[!is.na(MTF_2022$grade),]
table(MTF_2022$grade, useNA="always")   #grade

```

```{r, echo = T}
# 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")
```

## Calculate survey metrics

```{r, echo = T}
# 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
)
```

### Overall survey metrics

```{r, echo = T}
# 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
```

### Grouped by sex

```{r, echo = T}
# 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
```

### Grouped by grade and sex

```{r,echo=T}
# 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
```

## Process Data

```{r,echo=T}
# 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)
```

# NYTS
**National Youth Tobacco Survey (NYTS)**

## Prepare data
```{r, echo = T}
# 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")

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

# Age: 12-14, 15-17, 18+
table(NYTS_2023$age_cat, useNA = "always")

# Sex: Male, Female
table(NYTS_2023$um_sex_lab, useNA = "always")

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

```{r, echo = T}
# 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")
```

## Calculate survey metrics

```{r, echo = T}
# 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
)
```

### Overall survey metrics

```{r, echo = T}
# 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
```

### Grouped by age
```{r, echo = T}
# 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
```

### Grouped by age and sex
```{r,echo=T}
# 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
```

## Process Data

```{r,echo=T}
# 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)
```

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

## Prepare data

```{r, echo = T}
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")

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

# Age: 1=18-24, 2=25-34, 3=35-44, 4=45-54, 5=55-64, 6=65+
table(TUSCPS_2022$Age, useNA = "always")

# Sex:1=Male, 2=Female
table(TUSCPS_2022$Sex, useNA = "always")

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

```{r, echo = T}
# 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")
```

## Calculate survey metrics

```{r, echo = T}
# 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
  )
```

### Overall survey metrics

```{r, echo = T}
# 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
```

### Grouped by age

```{r, echo = T}
# 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
```

### Grouped by age and sex

```{r,echo=T}
# 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
```

## Process Data

```{r,echo=T}
# 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)
```


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

## Prepare data

```{r, echo = T}
# 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")

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

# Age: ~12, 13, 14, 15, 16, 17, 18+, Missing
table(YRBS_2021$Q1, useNA = "always")

# Sex: Male, Female, Missing
table(YRBS_2021$sex, useNA = "always")

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

YRBS_2021 <- YRBS_2021[YRBS_2021$sex != "Missing",]  # sex
YRBS_2021 <- YRBS_2021[!is.na(YRBS_2021$sex),]
table(YRBS_2021$sex, useNA="always")
```

```{r, echo = T}
# 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")
```

## Calculate survey metrics

```{r, echo = T}
# 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
)
```


### Overall survey metrics

```{r, echo = T}
# 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
```

### Grouped by age

```{r, echo = T}
# 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
```

### Grouped by age and sex

```{r,echo=T}
# 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
```

## Process Data

```{r,echo=T}
# 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)
```
