Load Packages

library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.1     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.5.0     ✔ tibble    3.2.1
## ✔ lubridate 1.9.2     ✔ tidyr     1.3.0
## ✔ purrr     1.0.1     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(janitor)
## 
## 载入程辑包：'janitor'
## 
## The following objects are masked from 'package:stats':
## 
##     chisq.test, fisher.test
library(rstudioapi)
library(docstring)
## 
## 载入程辑包：'docstring'
## 
## The following object is masked from 'package:utils':
## 
##     ?
library(dbplyr)
## 
## 载入程辑包：'dbplyr'
## 
## The following objects are masked from 'package:dplyr':
## 
##     ident, sql
library(purrr)
library(stringr)
library(reshape2)
## 
## 载入程辑包：'reshape2'
## 
## The following object is masked from 'package:tidyr':
## 
##     smiths
library(data.table)
## 
## 载入程辑包：'data.table'
## 
## The following objects are masked from 'package:reshape2':
## 
##     dcast, melt
## 
## The following objects are masked from 'package:lubridate':
## 
##     hour, isoweek, mday, minute, month, quarter, second, wday, week,
##     yday, year
## 
## The following objects are masked from 'package:dplyr':
## 
##     between, first, last
## 
## The following object is masked from 'package:purrr':
## 
##     transpose
library(knitr)
library(kableExtra)
## 
## 载入程辑包：'kableExtra'
## 
## The following object is masked from 'package:dplyr':
## 
##     group_rows

Part 1: Wrangle the data

For Part 1, the primary challenge was reshaping the dataset from a wide to a long format. The initial instructions lacked sufficient details about the required output format, so I reached out to Kevin for additional information. Based on his input, I selected the following outcomes from the original cleaned_quit_and_co.csv: “y1_tppd” (year 1 tobacco products per day), “y2_tppd”, “y3_tppd”, “y1_CO_recoded”, “y2_CO_recoded”, and “y3_CO_recoded”. I utilized dbplyr for dataset manipulation due to its SQL-like logic. Additionally, I converted the final output variables “tppd” and “CO” to numeric data types in preparation for Part 2.

The next step involved merging the “baseline.csv” with my newly created dataset “dat_split1”. By renaming the “study_id” column in the baseline dataset to “studyid”, I was able to complete this part with a simple merge function.

# 1.2) Set Working Directory to local version (set to data folder)

wdpath <- selectDirectory(caption = "Select Working Directory", label = "Set WD")
setwd(wdpath)

# 1.3) Load Dataset

dat_baseline <- read.csv("baseline.csv", header = TRUE)
dat_cleaned <- read.csv("cleaned_quit_and_co.csv", header = TRUE)

# Check the structure of the cleaned dataset
str(dat_cleaned)
## 'data.frame':    1166 obs. of  90 variables:
##  $ studyid                  : int  8053 7114 9158 8813 7297 8593 8468 7043 8698 8815 ...
##  $ TAU                      : int  1 NA NA 1 1 NA NA NA 1 1 ...
##  $ CHW_rand                 : int  0 1 1 0 0 1 1 1 0 0 ...
##  $ Arm                      : chr  "TAU (Treatment as usual)" "CHW only (Community health worker)" "AD + CHW" "TAU (Treatment as usual)" ...
##  $ Final_visit              : chr  "Baseline" "Baseline" "Year 3" "Year 3" ...
##  $ PCP.B.NPI                : int  730 757 401 489 229 730 652 295 664 759 ...
##  $ PCP.B.Name               : int  730 757 401 489 229 730 NA 295 664 759 ...
##  $ PCP.B.Clinic             : int  268 170 26 165 63 268 16 260 190 190 ...
##  $ PCP.B.Clinic_detailed    : chr  "No" "Yes" "No" "No" ...
##  $ PCP.B.Detailed           : logi  FALSE FALSE FALSE FALSE FALSE FALSE ...
##  $ PCP.B.Times_detailed     : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ PCP.Y2.NPI               : int  NA NA 401 755 229 541 NA 636 664 296 ...
##  $ PCP.Y2.Name              : int  NA NA 401 755 229 541 NA 636 664 296 ...
##  $ PCP.Y2.Clinic            : int  16 16 26 173 63 77 16 187 190 109 ...
##  $ PCP.Y2.Clinic_detailed   : chr  " " " " "No" "No" ...
##  $ PCP.Y2.Detailed          : logi  FALSE FALSE FALSE FALSE FALSE FALSE ...
##  $ PCP.Y2.Times_detailed    : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ PCP.Y3.NPI               : int  NA NA 401 755 220 541 NA 636 NA 296 ...
##  $ PCP.Y3.Name              : int  NA NA 401 755 220 541 NA 636 NA 296 ...
##  $ PCP.Y3.Clinic            : int  16 16 26 173 63 211 16 187 16 109 ...
##  $ PCP.Y3.Clinic_detailed   : chr  " " " " "No" "No" ...
##  $ PCP.Y3.Detailed          : chr  "FALSE" "FALSE" "FALSE" "FALSE" ...
##  $ PCP.Y3.Times_detailed    : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ PCP.Changed_from_Y2_to_Y3: chr  "No year 3" "No year 3" "No" "No" ...
##  $ PCP.Any_PCP_detailed     : chr  "FALSE" "FALSE" "FALSE" "FALSE" ...
##  $ Date.Baseline            : chr  "3/29/2017" "11/9/2016" "1/31/2017" "5/11/2017" ...
##  $ Date.Y2                  : chr  " " " " "4/15/2019" "10/2/2018" ...
##  $ Date.Y3                  : chr  " " " " "11/25/2019" "9/18/2019" ...
##  $ SCM.Y2.Any_prescribed    : chr  " " " " "No" "No" ...
##  $ SCM.Y3.Any_prescribed    : chr  " " " " "No" "No" ...
##  $ CHW.Consented            : chr  " " " " " " "Yes" ...
##  $ CHW.More_than_1_visit    : chr  " " " " " " "Yes" ...
##  $ CHW.Number_of_visits     : int  NA NA NA 12 NA NA NA 68 NA NA ...
##  $ SSS.CHR.Randomization    : chr  " " "CHW only (Community health worker)" "AD + CHW" " " ...
##  $ CHW.INT.Visits           : int  NA NA NA 12 NA NA NA 64 NA NA ...
##  $ CHW.INT.CBT_visits       : int  NA NA NA 0 NA NA NA 10 NA NA ...
##  $ rx_FINAL                 : chr  " " " " "CHW" " " ...
##  $ group_FINAL              : chr  " " " " "AD" " " ...
##  $ engagementqual           : chr  " " " " "Minimally engaged with CHW" " " ...
##  $ Engaged_BayCoveOnly      : int  NA NA 1 NA NA NA 0 NA NA NA ...
##  $ smallclinicrandomizarion : int  NA 1 NA NA NA NA NA NA NA NA ...
##  $ CBT_GROUPS_1             : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ CBT_GROUPS_2TO5          : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ CBT_GROUPS_6TO10         : int  NA NA NA NA NA NA NA 1 NA NA ...
##  $ CBT_GROUPS_11TO20        : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ CBT_GROUPS_21PLUS        : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ CBT_GROUPS_2ORMORE       : int  NA NA NA NA NA NA NA 1 NA NA ...
##  $ PCPEVERDETAILED          : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ GE1_CBTVISIT             : int  NA NA NA 0 NA NA NA 1 NA NA ...
##  $ HousingType              : chr  " " "Unknown" "24 Hour Group Home" "Independent" ...
##  $ VF_Active_OCT_2019       : chr  " " " " "Active" "Active" ...
##  $ GroupSharedLIving        : int  0 0 1 0 0 0 1 0 0 1 ...
##  $ ANYTAU                   : int  1 NA NA 1 1 NA NA NA 1 1 ...
##  $ y1_date                  : chr  "3/29/17 9:43" "11/9/16 10:49" "1/31/17 9:38" "5/11/17 12:21" ...
##  $ y1_tppd                  : int  4 16 36 18 8 5 40 20 23 16 ...
##  $ y1_tobaccotype_notes     : chr  "Cigarettes - Newports" "Cigarettes" "Cigarettes" "Cigarettes (Pall mall L&M)" ...
##  $ y1_cigs                  : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ y1_cigars                : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ y1_rollies               : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ y1_ecigsvape             : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ y1_morethanone           : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ y1_tobaccotypes          : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ y1_CO_recoded            : int  2 104 20 21 6 17 19 30 47 20 ...
##  $ Y1_CO_NOTES              : chr  " " " " " " " " ...
##  $ y1_CO_time               : chr  "3/29/17 9:49" "1/4/17 12:17" "1/31/17 9:47" "5/11/17 12:31" ...
##  $ y1_lastsmoke             : chr  "3/28/17 3:00" "1/4/17 12:10" "1/31/17 8:44" "5/11/17 11:50" ...
##  $ Y2_QUIT                  : int  NA NA 0 0 0 1 0 0 0 0 ...
##  $ y2_tppd                  : num  NA NA 27 12 0.5 0 40 6 4 20 ...
##  $ y2_CO_recoded            : chr  NA NA "43" "23" ...
##  $ Y2_CO_NOTES              : chr  " " " " " " " " ...
##  $ y2_co_time               : chr  NA NA "4/15/19 15:55" "10/2/18 9:57" ...
##  $ y2_lastsmoke             : chr  NA NA "4/15/19 15:25" "10/2/18 8:30" ...
##  $ y2_date                  : chr  NA NA "1/31/17" "5/11/17" ...
##  $ Y3_QUIT                  : int  NA NA 0 0 0 0 NA 0 NA 0 ...
##  $ y3_tppd                  : chr  NA NA "30" "20" ...
##  $ Y3_CO_NOTES              : chr  " " " " " " " " ...
##  $ y3_CO_recoded            : int  NA NA 51 38 NA 2 NA 21 NA 45 ...
##  $ y3_CO_time               : chr  NA NA "43794.66667" "9/18/19 15:29" ...
##  $ y3_lastsmoke             : chr  NA NA "11/25/19 16:17" "9/18/19 15:40" ...
##  $ y3_lastsmoke_updated     : chr  NA NA " " " " ...
##  $ y3_date                  : chr  NA NA "11/25/19 16:10" "9/18/19 15:29" ...
##  $ y3_number_product_types  : int  NA NA 1 1 1 1 NA 1 NA 1 ...
##  $ y3_concurrent            : int  NA NA 999 999 999 999 NA 999 NA 999 ...
##  $ y3_cigs                  : chr  NA NA "30" " " ...
##  $ y3_minicigs              : chr  NA NA " " " " ...
##  $ y3_rollies               : chr  NA NA " " "20" ...
##  $ y3_pipe                  : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ y3_ecigsvape             : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ y3_other                 : chr  NA NA " " " " ...
##  $ filter_.                 : int  0 0 0 0 0 0 0 0 0 0 ...
## Warning: Expected 2 pieces. Additional pieces discarded in 3498 rows [3499, 3500, 3501,
## 3502, 3503, 3504, 3505, 3506, 3507, 3508, 3509, 3510, 3511, 3512, 3513, 3514,
## 3515, 3516, 3517, 3518, ...].
## Warning: 强制改变过程中产生了NA

## Warning: 强制改变过程中产生了NA

Part 2: Summarize the data

This task presented two main challenges. First, it required splitting the new dataset “dat_merged” into a sample group randomly. I approached this by selecting 40% of the distinct study IDs and including all rows containing these IDs in the original dataset, resulting in the “dat_final” dataset ready for Part 3. Alternatively, one could perform the group split during Part 1 with the wide format and then convert it to the long format. Given that this task was in Part 2, I used the study ID selection strategy.

Second, it was necessary to include both gender and race information in one demographic summary. I converted all race codes to their respective labels as per the codebook, then split “dat_final” into “male” and “female” subsets and summarized each. This provided three types of summary tables: male, female, and overall. I then counted the number of each race within these subsets, creating a demographic table with the help of Copilot (a local AI assistant). Finally, I used the kable function to enhance the table’s presentation.

# 3.1) Summarize statistics

dat_grouped <- dat_merged %>%
  select("studyid", "Arm", "time", "age_1", "sex_1", "race_1", "tppd", "CO") %>%
  mutate(time = recode(time, "y1" = "Year 1", "y2" = "Year 2", "y3" = "Year 3"))

# 3.2) Split dataset to 2 randomization group

# Randomly select 40% of unique studyid
set.seed(123)

selected_studyid <- dat_grouped %>%
  distinct(studyid) %>%
  sample_frac(0.4) %>%
  pull(studyid)

# Filter the dataset to include all records for the selected studyid
dat_final <- dat_grouped %>%
  filter(studyid %in% selected_studyid)


# 3.3) Split dataset by gender

male <- dat_final %>%
  filter(sex_1 == 1)

female <- dat_final %>%
  filter(sex_1 == 2)


# 3.4) Count race

# Count the number of each race in the overall dataset
race_counts_overall <- table(dat_final$race_1)

# Count the number of each race in the male subset
race_counts_male <- table(male$race_1)

# Count the number of each race in the female subset
race_counts_female <- table(female$race_1)

Create the demographic table manually. This approach is more flexible to edit the demograpic table and saves time with Copilot’s assistance, even though some packages can automate it.

# 3.5) Create demographic summary

# Manually typed summary table
race_labels <- c("White", "Haitian, Black, or African American", "Asian",
                 "Hawaiian or other Pacific Islander", "American Indian/Alaska Native",
                 "Other", "Multi-race")

# Provided counts
race_counts_overall <- c(621, 507, 54, 3, 12, 102, 99)
race_counts_male <- c(405, 369, 45, 0, 6, 72, 72)
race_counts_female <- c(216, 138, 9, 3, 6, 30, 27)

# Create the demographic summary data frame
demographic_summary <- data.frame(
  Metric = c("Count", "Age Mean", "Age SD", "TPPD Mean", "TPPD SD", "Expired CO Mean", "Expired CO SD", "NA Count of TPPD", "NA Count of Expired CO",
             paste("Race", race_labels, "Count", sep = " - ")),
  Male = c(
    nrow(male),                          # Count of males
    mean(male$age_1, na.rm = TRUE),      # Mean age of males
    sd(male$age_1, na.rm = TRUE),        # Standard deviation of age of males
    mean(male$tppd, na.rm = TRUE),       # Mean TPPD for males
    sd(male$tppd, na.rm = TRUE),         # Standard deviation of TPPD for males
    mean(male$CO, na.rm = TRUE),         # Mean Expired CO for males
    sd(male$CO, na.rm = TRUE),           # Standard deviation of Expired CO for males
    sum(is.na(male$tppd)),               # Number of NA in TPPD for males
    sum(is.na(male$CO)),                 # Number of NA in Expired CO for males
    race_counts_male                     # Race counts for males
  ),
  Female = c(
    nrow(female),                        # Count of females
    mean(female$age_1, na.rm = TRUE),    # Mean age of females
    sd(female$age_1, na.rm = TRUE),      # Standard deviation of age of females
    mean(female$tppd, na.rm = TRUE),     # Mean TPPD for females
    sd(female$tppd, na.rm = TRUE),       # Standard deviation of TPPD for females
    mean(female$CO, na.rm = TRUE),       # Mean Expired CO for females
    sd(female$CO, na.rm = TRUE),         # Standard deviation of Expired CO for females
    sum(is.na(female$tppd)),             # Number of NA in TPPD for females
    sum(is.na(female$CO)),               # Number of NA in Expired CO for females
    race_counts_female                   # Race counts for females
  ),
  Overall = c(
    nrow(dat_final),                           # Count of overall dataset
    mean(dat_final$age_1, na.rm = TRUE),       # Mean age of overall dataset
    sd(dat_final$age_1, na.rm = TRUE),         # Standard deviation of age of overall dataset
    mean(dat_final$tppd, na.rm = TRUE),        # Mean TPPD for overall dataset
    sd(dat_final$tppd, na.rm = TRUE),          # Standard deviation of TPPD for overall dataset
    mean(dat_final$CO, na.rm = TRUE),          # Mean Expired CO for overall dataset
    sd(dat_final$CO, na.rm = TRUE),            # Standard deviation of Expired CO for overall dataset
    sum(is.na(dat_final$tppd)),                # Number of NA in TPPD for overall dataset
    sum(is.na(dat_final$CO)),                  # Number of NA in Expired CO for overall dataset
    race_counts_overall                        # Race counts for overall dataset
  )
)

# Print the demographic summary table
kable(demographic_summary[-1, ], format = "html", col.names = c("", "Male", "Female", "Overall")) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = F) %>%
  column_spec(1, bold = TRUE) %>%
  add_header_above(c("Demographic Summary" = 5))
Demographic Summary
Male Female Overall
2 Age Mean 47.37461 47.69231 47.47210
3 Age SD 13.23356 13.72386 13.38180
4 TPPD Mean 14.48114 14.12791 14.37318
5 TPPD SD 11.32228 10.73410 11.14210
6 Expired CO Mean 25.59192 22.00000 24.49903
7 Expired CO SD 21.74229 16.27637 20.29568
8 NA Count of TPPD 208.00000 94.00000 302.00000
9 NA Count of Expired CO 251.00000 115.00000 366.00000
10 Race - White - Count 405.00000 216.00000 621.00000
11 Race - Haitian, Black, or African American - Count 369.00000 138.00000 507.00000
12 Race - Asian - Count 45.00000 9.00000 54.00000
13 Race - Hawaiian or other Pacific Islander - Count 0.00000 3.00000 3.00000
14 Race - American Indian/Alaska Native - Count 6.00000 6.00000 12.00000
15 Race - Other - Count 72.00000 30.00000 102.00000
16 Race - Multi-race - Count 72.00000 27.00000 99.00000

Part 3: Graph the data

Part 3 was somewhat challenging. I had to decide which outcome to use for the sole plot submission. Initially, I considered using “tppd” (tobacco products per day) as the outcome. However, after Kevin explained the significance of “CO” (expired CO), I realized it would be a more objective outcome compared to “tppd” (which seemed self-reported).

# 4.1) Plot Expired CO Distribution

plot <- ggplot(dat_final, aes(x = time, y = CO, fill = Arm)) +
  geom_boxplot() +
  labs(title = "Expired CO Distribution by Experimental Arm and Time",
       x = "Time",
       y = "Expired CO",
       fill = "Experimental Arm")

# Print the plot
print(plot)
## Warning: Removed 366 rows containing non-finite outside the scale range
## (`stat_boxplot()`).

The boxplot showing “Expired CO” values over time for different experimental groups is good for seeing the overall effects of the interventions. However, the “AD + CHW” and “AD only” groups don’t show much change over time. To better understand these small differences, we could use violin plots or other types of plots. These might show more details about the data and help us find important patterns that the boxplot doesn’t show.