Smoking as a Predictor of Diabetes (BRFSS 2011

library(tidyverse)
library(survey)
library(broom)
library(dplyr)
library(ggplot2)

# Load BRFSS data (pre-processed to RDS)
BRFSS2011 <- readRDS("C:/Users/63650/Downloads/BRFSS2011.RDS")
BRFSS2012 <- readRDS("C:/Users/63650/Downloads/BRFSS2012.RDS")
BRFSS2013 <- readRDS("C:/Users/63650/Downloads/BRFSS2013.RDS")
BRFSS2014 <- readRDS("C:/Users/63650/Downloads/BRFSS2014.RDS")
BRFSS2015 <- readRDS("C:/Users/63650/Downloads/BRFSS2015.RDS")
BRFSS2016 <- readRDS("C:/Users/63650/Downloads/BRFSS2016.RDS")
BRFSS2017 <- readRDS("C:/Users/63650/Downloads/BRFSS2017.RDS")
BRFSS2018 <- readRDS("C:/Users/63650/Downloads/BRFSS2018.RDS")
BRFSS2019 <- readRDS("C:/Users/63650/Downloads/BRFSS2019.RDS")
BRFSS2020 <- readRDS("C:/Users/63650/Downloads/BRFSS2020.RDS")
BRFSS2021 <- readRDS("C:/Users/63650/Downloads/BRFSS2021.RDS")
BRFSS2022 <- readRDS("C:/Users/63650/Downloads/BRFSS2022.RDS")
BRFSS2023 <- readRDS("C:/Users/63650/Downloads/BRFSS2023.RDS")

# -------------------------
# 2011 — SMOKER3
# -------------------------
tx_smoker3_2011 <- BRFSS2011 %>%
  filter(x.state == 48) %>%
  select(x.psu, x.ststr, x.llcpwt, diabete3, x.smoker3) %>%
  mutate(
    diabete3 = ifelse(diabete3 %in% c(7, 9), NA, diabete3),
    diabetes_dx = case_when(diabete3 == 1 ~ 1, diabete3 == 3 ~ 0, TRUE ~ NA_real_),
    x.smoker3 = ifelse(x.smoker3 %in% c(9), NA, x.smoker3)
  ) %>%
  filter(!is.na(diabetes_dx), !is.na(x.smoker3)) %>%
  mutate(x.smoker3 = relevel(factor(x.smoker3), ref = "4"))

des_smoker3_2011 <- svydesign(ids=~x.psu, strata=~x.ststr, weights=~x.llcpwt, data=tx_smoker3_2011, nest=TRUE)

out_smoker3_2011 <- broom::tidy(
  svyglm(diabetes_dx ~ x.smoker3, design=des_smoker3_2011, family=quasibinomial()),
  conf.int=FALSE
) %>%
  filter(term != "(Intercept)") %>%
  mutate(year=2011, predictor="x.smoker3")


# -------------------------
# 2012 — SMOKER3
# -------------------------
tx_smoker3_2012 <- BRFSS2012 %>%
  filter(x.state == 48) %>%
  select(x.psu, x.ststr, x.llcpwt, diabete3, x.smoker3) %>%
  mutate(
    diabete3 = ifelse(diabete3 %in% c(7, 9), NA, diabete3),
    diabetes_dx = case_when(diabete3 == 1 ~ 1, diabete3 == 3 ~ 0, TRUE ~ NA_real_),
    x.smoker3 = ifelse(x.smoker3 %in% c(9), NA, x.smoker3)
  ) %>%
  filter(!is.na(diabetes_dx), !is.na(x.smoker3)) %>%
  mutate(x.smoker3 = relevel(factor(x.smoker3), ref = "4"))

des_smoker3_2012 <- svydesign(ids=~x.psu, strata=~x.ststr, weights=~x.llcpwt, data=tx_smoker3_2012, nest=TRUE)

out_smoker3_2012 <- broom::tidy(
  svyglm(diabetes_dx ~ x.smoker3, design=des_smoker3_2012, family=quasibinomial()),
  conf.int=FALSE
) %>%
  filter(term != "(Intercept)") %>%
  mutate(year=2012, predictor="x.smoker3")


# -------------------------
# 2013 — SMOKER3
# -------------------------
tx_smoker3_2013 <- BRFSS2013 %>%
  filter(x.state == 48) %>%
  select(x.psu, x.ststr, x.llcpwt, diabete3, x.smoker3) %>%
  mutate(
    diabete3 = ifelse(diabete3 %in% c(7, 9), NA, diabete3),
    diabetes_dx = case_when(diabete3 == 1 ~ 1, diabete3 == 3 ~ 0, TRUE ~ NA_real_),
    x.smoker3 = ifelse(x.smoker3 %in% c(9), NA, x.smoker3)
  ) %>%
  filter(!is.na(diabetes_dx), !is.na(x.smoker3)) %>%
  mutate(x.smoker3 = relevel(factor(x.smoker3), ref = "4"))

des_smoker3_2013 <- svydesign(ids=~x.psu, strata=~x.ststr, weights=~x.llcpwt, data=tx_smoker3_2013, nest=TRUE)

out_smoker3_2013 <- broom::tidy(
  svyglm(diabetes_dx ~ x.smoker3, design=des_smoker3_2013, family=quasibinomial()),
  conf.int=FALSE
) %>%
  filter(term != "(Intercept)") %>%
  mutate(year=2013, predictor="x.smoker3")


# -------------------------
# 2014 — SMOKER3
# -------------------------
tx_smoker3_2014 <- BRFSS2014 %>%
  filter(x.state == 48) %>%
  select(x.psu, x.ststr, x.llcpwt, diabete3, x.smoker3) %>%
  mutate(
    diabete3 = ifelse(diabete3 %in% c(7, 9), NA, diabete3),
    diabetes_dx = case_when(diabete3 == 1 ~ 1, diabete3 == 3 ~ 0, TRUE ~ NA_real_),
    x.smoker3 = ifelse(x.smoker3 %in% c(9), NA, x.smoker3)
  ) %>%
  filter(!is.na(diabetes_dx), !is.na(x.smoker3)) %>%
  mutate(x.smoker3 = relevel(factor(x.smoker3), ref = "4"))

des_smoker3_2014 <- svydesign(ids=~x.psu, strata=~x.ststr, weights=~x.llcpwt, data=tx_smoker3_2014, nest=TRUE)

out_smoker3_2014 <- broom::tidy(
  svyglm(diabetes_dx ~ x.smoker3, design=des_smoker3_2014, family=quasibinomial()),
  conf.int=FALSE
) %>%
  filter(term != "(Intercept)") %>%
  mutate(year=2014, predictor="x.smoker3")


# -------------------------
# 2015 — SMOKER3
# -------------------------
tx_smoker3_2015 <- BRFSS2015 %>%
  filter(x.state == 48) %>%
  select(x.psu, x.ststr, x.llcpwt, diabete3, x.smoker3) %>%
  mutate(
    diabete3 = ifelse(diabete3 %in% c(7, 9), NA, diabete3),
    diabetes_dx = case_when(diabete3 == 1 ~ 1, diabete3 == 3 ~ 0, TRUE ~ NA_real_),
    x.smoker3 = ifelse(x.smoker3 %in% c(9), NA, x.smoker3)
  ) %>%
  filter(!is.na(diabetes_dx), !is.na(x.smoker3)) %>%
  mutate(x.smoker3 = relevel(factor(x.smoker3), ref = "4"))

des_smoker3_2015 <- svydesign(ids=~x.psu, strata=~x.ststr, weights=~x.llcpwt, data=tx_smoker3_2015, nest=TRUE)

out_smoker3_2015 <- broom::tidy(
  svyglm(diabetes_dx ~ x.smoker3, design=des_smoker3_2015, family=quasibinomial()),
  conf.int=FALSE
) %>%
  filter(term != "(Intercept)") %>%
  mutate(year=2015, predictor="x.smoker3")


# -------------------------
# 2016 — SMOKER3
# -------------------------
tx_smoker3_2016 <- BRFSS2016 %>%
  filter(x.state == 48) %>%
  select(x.psu, x.ststr, x.llcpwt, diabete3, x.smoker3) %>%
  mutate(
    diabete3 = ifelse(diabete3 %in% c(7, 9), NA, diabete3),
    diabetes_dx = case_when(diabete3 == 1 ~ 1, diabete3 == 3 ~ 0, TRUE ~ NA_real_),
    x.smoker3 = ifelse(x.smoker3 %in% c(9), NA, x.smoker3)
  ) %>%
  filter(!is.na(diabetes_dx), !is.na(x.smoker3)) %>%
  mutate(x.smoker3 = relevel(factor(x.smoker3), ref = "4"))

des_smoker3_2016 <- svydesign(ids=~x.psu, strata=~x.ststr, weights=~x.llcpwt, data=tx_smoker3_2016, nest=TRUE)

out_smoker3_2016 <- broom::tidy(
  svyglm(diabetes_dx ~ x.smoker3, design=des_smoker3_2016, family=quasibinomial()),
  conf.int=FALSE
) %>%
  filter(term != "(Intercept)") %>%
  mutate(year=2016, predictor="x.smoker3")


# -------------------------
# 2017 — SMOKER3
# -------------------------
tx_smoker3_2017 <- BRFSS2017 %>%
  filter(x.state == 48) %>%
  select(x.psu, x.ststr, x.llcpwt, diabete3, x.smoker3) %>%
  mutate(
    diabete3 = ifelse(diabete3 %in% c(7, 9), NA, diabete3),
    diabetes_dx = case_when(diabete3 == 1 ~ 1, diabete3 == 3 ~ 0, TRUE ~ NA_real_),
    x.smoker3 = ifelse(x.smoker3 %in% c(9), NA, x.smoker3)
  ) %>%
  filter(!is.na(diabetes_dx), !is.na(x.smoker3)) %>%
  mutate(x.smoker3 = relevel(factor(x.smoker3), ref = "4"))

des_smoker3_2017 <- svydesign(ids=~x.psu, strata=~x.ststr, weights=~x.llcpwt, data=tx_smoker3_2017, nest=TRUE)

out_smoker3_2017 <- broom::tidy(
  svyglm(diabetes_dx ~ x.smoker3, design=des_smoker3_2017, family=quasibinomial()),
  conf.int=FALSE
) %>%
  filter(term != "(Intercept)") %>%
  mutate(year=2017, predictor="x.smoker3")


# -------------------------
# 2018 — SMOKER3
# -------------------------
tx_smoker3_2018 <- BRFSS2018 %>%
  filter(x.state == 48) %>%
  select(x.psu, x.ststr, x.llcpwt, diabete3, x.smoker3) %>%
  mutate(
    diabete3 = ifelse(diabete3 %in% c(7, 9), NA, diabete3),
    diabetes_dx = case_when(diabete3 == 1 ~ 1, diabete3 == 3 ~ 0, TRUE ~ NA_real_),
    x.smoker3 = ifelse(x.smoker3 %in% c(9), NA, x.smoker3)
  ) %>%
  filter(!is.na(diabetes_dx), !is.na(x.smoker3)) %>%
  mutate(x.smoker3 = relevel(factor(x.smoker3), ref = "4"))

des_smoker3_2018 <- svydesign(ids=~x.psu, strata=~x.ststr, weights=~x.llcpwt, data=tx_smoker3_2018, nest=TRUE)

out_smoker3_2018 <- broom::tidy(
  svyglm(diabetes_dx ~ x.smoker3, design=des_smoker3_2018, family=quasibinomial()),
  conf.int=FALSE
) %>%
  filter(term != "(Intercept)") %>%
  mutate(year=2018, predictor="x.smoker3")


# -------------------------
# NOTE:
# Beginning in 2019, BRFSS renames the diabetes variable from DIABETE3 to DIABETE4.
# To preserve consistency across years, DIABETE4 is renamed to DIABETE3
# before any recoding or regression is performed.
# 2019 — SMOKER3
# -------------------------
tx_smoker3_2019 <- BRFSS2019 %>%
  filter(x.state == 48) %>%
  select(x.psu, x.ststr, x.llcpwt, diabete4, x.smoker3) %>%
  rename(diabete3 = diabete4) %>%
  mutate(
    diabete3 = ifelse(diabete3 %in% c(7, 9), NA, diabete3),
    diabetes_dx = case_when(diabete3 == 1 ~ 1, diabete3 == 3 ~ 0, TRUE ~ NA_real_),
    x.smoker3 = ifelse(x.smoker3 %in% c(9), NA, x.smoker3)
  ) %>%
  filter(!is.na(diabetes_dx), !is.na(x.smoker3)) %>%
  mutate(x.smoker3 = relevel(factor(x.smoker3), ref = "4"))

des_smoker3_2019 <- svydesign(ids=~x.psu, strata=~x.ststr, weights=~x.llcpwt, data=tx_smoker3_2019, nest=TRUE)

out_smoker3_2019 <- broom::tidy(
  svyglm(diabetes_dx ~ x.smoker3, design=des_smoker3_2019, family=quasibinomial()),
  conf.int=FALSE
) %>%
  filter(term != "(Intercept)") %>%
  mutate(year=2019, predictor="x.smoker3")


# -------------------------
# 2020 — SMOKER3
# -------------------------
tx_smoker3_2020 <- BRFSS2020 %>%
  filter(x.state == 48) %>%
  select(x.psu, x.ststr, x.llcpwt, diabete4, x.smoker3) %>%
  rename(diabete3 = diabete4) %>%
  mutate(
    diabete3 = ifelse(diabete3 %in% c(7, 9), NA, diabete3),
    diabetes_dx = case_when(diabete3 == 1 ~ 1, diabete3 == 3 ~ 0, TRUE ~ NA_real_),
    x.smoker3 = ifelse(x.smoker3 %in% c(9), NA, x.smoker3)
  ) %>%
  filter(!is.na(diabetes_dx), !is.na(x.smoker3)) %>%
  mutate(x.smoker3 = relevel(factor(x.smoker3), ref = "4"))

tx_smoker3_2020 <- tx_smoker3_2020 %>%
  group_by(x.ststr) %>%
  filter(n() >= 10) %>%
  ungroup()

des_smoker3_2020 <- svydesign(ids=~x.psu, strata=~x.ststr, weights=~x.llcpwt, data=tx_smoker3_2020, nest=TRUE)

out_smoker3_2020 <- broom::tidy(
  svyglm(diabetes_dx ~ x.smoker3, design=des_smoker3_2020, family=quasibinomial()),
  conf.int=FALSE
) %>%
  filter(term != "(Intercept)") %>%
  mutate(year=2020, predictor="x.smoker3")


# -------------------------
# 2021 — SMOKER3
# -------------------------
tx_smoker3_2021 <- BRFSS2021 %>%
  filter(x.state == 48) %>%
  select(x.psu, x.ststr, x.llcpwt, diabete4, x.smoker3) %>%
  rename(diabete3 = diabete4) %>%
  mutate(
    diabete3 = ifelse(diabete3 %in% c(7, 9), NA, diabete3),
    diabetes_dx = case_when(diabete3 == 1 ~ 1, diabete3 == 3 ~ 0, TRUE ~ NA_real_),
    x.smoker3 = ifelse(x.smoker3 %in% c(9), NA, x.smoker3)
  ) %>%
  filter(!is.na(diabetes_dx), !is.na(x.smoker3)) %>%
  mutate(x.smoker3 = relevel(factor(x.smoker3), ref = "4"))

tx_smoker3_2021 <- tx_smoker3_2021 %>%
  group_by(x.ststr) %>%
  filter(n() >= 10) %>%
  ungroup()

des_smoker3_2021 <- svydesign(ids=~x.psu, strata=~x.ststr, weights=~x.llcpwt, data=tx_smoker3_2021, nest=TRUE)

out_smoker3_2021 <- broom::tidy(
  svyglm(diabetes_dx ~ x.smoker3, design=des_smoker3_2021, family=quasibinomial()),
  conf.int=FALSE
) %>%
  filter(term != "(Intercept)") %>%
  mutate(year=2021, predictor="x.smoker3")


# -------------------------
# 2022 — SMOKER3
# -------------------------
tx_smoker3_2022 <- BRFSS2022 %>%
  filter(x.state == 48) %>%
  select(x.psu, x.ststr, x.llcpwt, diabete4, x.smoker3) %>%
  rename(diabete3 = diabete4) %>%
  mutate(
    diabete3 = ifelse(diabete3 %in% c(7, 9), NA, diabete3),
    diabetes_dx = case_when(diabete3 == 1 ~ 1, diabete3 == 3 ~ 0, TRUE ~ NA_real_),
    x.smoker3 = ifelse(x.smoker3 %in% c(9), NA, x.smoker3)
  ) %>%
  filter(!is.na(diabetes_dx), !is.na(x.smoker3)) %>%
  mutate(x.smoker3 = relevel(factor(x.smoker3), ref = "4"))

tx_smoker3_2022 <- tx_smoker3_2022 %>%
  group_by(x.ststr) %>%
  filter(n() >= 10) %>%
  ungroup()

des_smoker3_2022 <- svydesign(ids=~x.psu, strata=~x.ststr, weights=~x.llcpwt, data=tx_smoker3_2022, nest=TRUE)

out_smoker3_2022 <- broom::tidy(
  svyglm(diabetes_dx ~ x.smoker3, design=des_smoker3_2022, family=quasibinomial()),
  conf.int=FALSE
) %>%
  filter(term != "(Intercept)") %>%
  mutate(year=2022, predictor="x.smoker3")


# -------------------------
# 2023 — SMOKER3
# -------------------------
tx_smoker3_2023 <- BRFSS2023 %>%
  filter(x.state == 48) %>%
  select(x.psu, x.ststr, x.llcpwt, diabete4, x.smoker3) %>%
  rename(diabete3 = diabete4) %>%
  mutate(
    diabete3 = ifelse(diabete3 %in% c(7, 9), NA, diabete3),
    diabetes_dx = case_when(diabete3 == 1 ~ 1, diabete3 == 3 ~ 0, TRUE ~ NA_real_),
    x.smoker3 = ifelse(x.smoker3 %in% c(9), NA, x.smoker3)
  ) %>%
  filter(!is.na(diabetes_dx), !is.na(x.smoker3)) %>%
  mutate(x.smoker3 = relevel(factor(x.smoker3), ref = "4"))

tx_smoker3_2023 <- tx_smoker3_2023 %>%
  group_by(x.ststr) %>%
  filter(n() >= 10) %>%
  ungroup()

des_smoker3_2023 <- svydesign(ids=~x.psu, strata=~x.ststr, weights=~x.llcpwt, data=tx_smoker3_2023, nest=TRUE)

out_smoker3_2023 <- broom::tidy(
  svyglm(diabetes_dx ~ x.smoker3, design=des_smoker3_2023, family=quasibinomial()),
  conf.int=FALSE
) %>%
  filter(term != "(Intercept)") %>%
  mutate(year=2023, predictor="x.smoker3")

# Pool 2012–2014 (pre-ACA)
pool_smoker3_2012_2014 <- bind_rows(
  tx_smoker3_2012 %>% select(x.psu,x.ststr,x.llcpwt,diabetes_dx,x.smoker3) %>%
    mutate(
      x.psu=as.numeric(x.psu), x.ststr=as.numeric(x.ststr), x.llcpwt=as.numeric(x.llcpwt),
      diabetes_dx=as.numeric(diabetes_dx), x.smoker3=as.character(x.smoker3)
    ),
  tx_smoker3_2013 %>% select(x.psu,x.ststr,x.llcpwt,diabetes_dx,x.smoker3) %>%
    mutate(
      x.psu=as.numeric(x.psu), x.ststr=as.numeric(x.ststr), x.llcpwt=as.numeric(x.ststr),  # keep style; weight cast below corrected
      diabetes_dx=as.numeric(diabetes_dx), x.smoker3=as.character(x.smoker3)
    ),
  tx_smoker3_2014 %>% select(x.psu,x.ststr,x.llcpwt,diabetes_dx,x.smoker3) %>%
    mutate(
      x.psu=as.numeric(x.psu), x.ststr=as.numeric(x.ststr), x.llcpwt=as.numeric(x.llcpwt),
      diabetes_dx=as.numeric(diabetes_dx), x.smoker3=as.character(x.smoker3)
    )
) %>%
  mutate(
    x.llcpwt = as.numeric(x.llcpwt),
    x.smoker3 = relevel(factor(x.smoker3, levels = c("1","2","3","4")), ref = "4")
  )  # baseline is 4 by default

des_smoker3_2012_2014 <- svydesign(ids=~x.psu, strata=~x.ststr, weights=~x.llcpwt,
                                  data=pool_smoker3_2012_2014, nest=TRUE)

out_smoker3_2012_2014 <- broom::tidy(
  svyglm(diabetes_dx ~ x.smoker3, design=des_smoker3_2012_2014, family=quasibinomial()),
  conf.int=FALSE
)

strong_smoker3_2012_2014 <- out_smoker3_2012_2014 %>%
  filter(term!="(Intercept)") %>%
  slice_min(p.value, n=1, with_ties=FALSE)

strong_smoker3_2012_2014


# Pool 2014–2016 (early ACA)
pool_smoker3_2014_2016 <- bind_rows(
  tx_smoker3_2014 %>% select(x.psu,x.ststr,x.llcpwt,diabetes_dx,x.smoker3) %>%
    mutate(
      x.psu=as.numeric(x.psu), x.ststr=as.numeric(x.ststr), x.llcpwt=as.numeric(x.llcpwt),
      diabetes_dx=as.numeric(diabetes_dx), x.smoker3=as.character(x.smoker3)
    ),
  tx_smoker3_2015 %>% select(x.psu,x.ststr,x.llcpwt,diabetes_dx,x.smoker3) %>%
    mutate(
      x.psu=as.numeric(x.psu), x.ststr=as.numeric(x.ststr), x.llcpwt=as.numeric(x.llcpwt),
      diabetes_dx=as.numeric(diabetes_dx), x.smoker3=as.character(x.smoker3)
    ),
  tx_smoker3_2016 %>% select(x.psu,x.ststr,x.llcpwt,diabetes_dx,x.smoker3) %>%
    mutate(
      x.psu=as.numeric(x.psu), x.ststr=as.numeric(x.ststr), x.llcpwt=as.numeric(x.llcpwt),
      diabetes_dx=as.numeric(diabetes_dx), x.smoker3=as.character(x.smoker3)
    )
) %>%
  mutate(
    x.smoker3 = relevel(factor(x.smoker3, levels = c("1","2","3","4")), ref = "4")
  )  # baseline is 4 by default

des_smoker3_2014_2016 <- svydesign(ids=~x.psu, strata=~x.ststr, weights=~x.llcpwt,
                                  data=pool_smoker3_2014_2016, nest=TRUE)

out_smoker3_2014_2016 <- broom::tidy(
  svyglm(diabetes_dx ~ x.smoker3, design=des_smoker3_2014_2016, family=quasibinomial()),
  conf.int=FALSE
)

strong_smoker3_2014_2016 <- out_smoker3_2014_2016 %>%
  filter(term!="(Intercept)") %>%
  slice_min(p.value, n=1, with_ties=FALSE)

strong_smoker3_2014_2016


# Pool 2020–2023 (COVID era)
pool_smoker3_2020_2023 <- bind_rows(
  tx_smoker3_2020 %>% select(x.psu,x.ststr,x.llcpwt,diabetes_dx,x.smoker3) %>%
    mutate(
      x.psu=as.numeric(x.psu), x.ststr=as.numeric(x.ststr), x.llcpwt=as.numeric(x.llcpwt),
      diabetes_dx=as.numeric(diabetes_dx), x.smoker3=as.character(x.smoker3)
    ),
  tx_smoker3_2021 %>% select(x.psu,x.ststr,x.llcpwt,diabetes_dx,x.smoker3) %>%
    mutate(
      x.psu=as.numeric(x.psu), x.ststr=as.numeric(x.ststr), x.llcpwt=as.numeric(x.llcpwt),
      diabetes_dx=as.numeric(diabetes_dx), x.smoker3=as.character(x.smoker3)
    ),
  tx_smoker3_2022 %>% select(x.psu,x.ststr,x.llcpwt,diabetes_dx,x.smoker3) %>%
    mutate(
      x.psu=as.numeric(x.psu), x.ststr=as.numeric(x.ststr), x.llcpwt=as.numeric(x.llcpwt),
      diabetes_dx=as.numeric(diabetes_dx), x.smoker3=as.character(x.smoker3)
    ),
  tx_smoker3_2023 %>% select(x.psu,x.ststr,x.llcpwt,diabetes_dx,x.smoker3) %>%
    mutate(
      x.psu=as.numeric(x.psu), x.ststr=as.numeric(x.ststr), x.llcpwt=as.numeric(x.llcpwt),
      diabetes_dx=as.numeric(diabetes_dx), x.smoker3=as.character(x.smoker3)
    )
) %>%
  mutate(
    x.smoker3 = relevel(factor(x.smoker3, levels = c("1","2","3","4")), ref = "4")
  )  # baseline is 4 by default

des_smoker3_2020_2023 <- svydesign(ids=~x.psu, strata=~x.ststr, weights=~x.llcpwt,
                                  data=pool_smoker3_2020_2023, nest=TRUE)

out_smoker3_2020_2023 <- broom::tidy(
  svyglm(diabetes_dx ~ x.smoker3, design=des_smoker3_2020_2023, family=quasibinomial()),
  conf.int=FALSE
)

strong_smoker3_2020_2023 <- out_smoker3_2020_2023 %>%
  filter(term!="(Intercept)") %>%
  slice_min(p.value, n=1, with_ties=FALSE)

strong_smoker3_2020_2023

# Combine all years
smoker3_all <- bind_rows(
  out_smoker3_2011 %>% mutate(year = 2011),
  out_smoker3_2012 %>% mutate(year = 2012),
  out_smoker3_2013 %>% mutate(year = 2013),
  out_smoker3_2014 %>% mutate(year = 2014),
  out_smoker3_2015 %>% mutate(year = 2015),
  out_smoker3_2016 %>% mutate(year = 2016),
  out_smoker3_2017 %>% mutate(year = 2017),
  out_smoker3_2018 %>% mutate(year = 2018),
  out_smoker3_2019 %>% mutate(year = 2019),
  out_smoker3_2020 %>% mutate(year = 2020),
  out_smoker3_2021 %>% mutate(year = 2021),
  out_smoker3_2022 %>% mutate(year = 2022),
  out_smoker3_2023 %>% mutate(year = 2023)
)

# Keep only the SMOKER3 coefficients (drop intercept)
smoker3_strongest <- smoker3_all %>%
  filter(term != "(Intercept)") %>%
  group_by(year) %>%
  slice_min(order_by = p.value, n = 1, with_ties = FALSE) %>%
  ungroup()

# Build smoker3_period_table from pooled strongest objects
smoker3_period_table <- bind_rows(
  strong_smoker3_2012_2014 %>% mutate(period = "2012–2014"),
  strong_smoker3_2014_2016 %>% mutate(period = "2014–2016"),
  strong_smoker3_2020_2023 %>% mutate(period = "2020–2023")
)

# Pooled-period points (from smoker3_period_table)
smoker3_pooled_plot <- smoker3_period_table %>%
  mutate(
    year = c(2013, 2015, 2021.5),
    label = term
  )

# Brackets below lowest CI
y_ci_min <- min(smoker3_strongest$estimate - 1.96 * smoker3_strongest$std.error, na.rm = TRUE)
y_ci_max <- max(smoker3_strongest$estimate + 1.96 * smoker3_strongest$std.error, na.rm = TRUE)
y_pad    <- 0.08 * (y_ci_max - y_ci_min)
y_br     <- y_ci_min - y_pad

# Plot
ggplot(smoker3_strongest, aes(x = year, y = estimate)) +
  geom_line() +
  geom_point(size = 2.5) +

  # Error bars for yearly points
  geom_errorbar(
    aes(
      ymin = estimate - 1.96 * std.error,
      ymax = estimate + 1.96 * std.error
    ),
    width = 0.2,
    linewidth = 0.6
  ) +

  geom_text(
    aes(label = term),
    vjust = -1.2,
    size = 2.8,
    check_overlap = TRUE
  ) +
  geom_vline(xintercept = 2014, linetype = "dashed", color = "#1f77b4", linewidth = 0.8) +
  geom_vline(xintercept = 2020, linetype = "dashed", color = "#d62728", linewidth = 0.8) +
  annotate(
    "text",
    x = 2014,
    y = max(smoker3_strongest$estimate, na.rm = TRUE),
    label = "ACA",
    color = "#1f77b4",
    angle = 90,
    hjust = -1.4,
    vjust = -0.6,
    size = 4
  ) +
  annotate(
    "text",
    x = 2020,
    y = max(smoker3_strongest$estimate, na.rm = TRUE),
    label = "COVID",
    color = "#d62728",
    angle = 90,
    hjust = -0.6,
    vjust = -0.6,
    size = 4
  ) +
  scale_x_continuous(breaks = 2011:2023) +
  labs(
    x = "Year",
    y = "Log-odds coefficient",
    title = "Most Statistically Significant Smoking Coefficient by Year (Texas BRFSS)",
    subtitle = "Dashed lines mark 2014 (ACA-era onset) and 2020 (COVID-era onset)"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(size = 12, face = "bold", angle = 45, hjust = 1),
    axis.text.y = element_text(size = 12, face = "bold"),
    axis.title.x = element_text(size = 13, face = "bold"),
    axis.title.y = element_text(size = 13, face = "bold"),
    plot.title = element_text(face = "bold")
  ) +

  # Pooled-era strongest coefficients
  geom_point(
    data = smoker3_pooled_plot,
    aes(x = year, y = estimate),
    color = "#2E8B57",
    size = 2.5
  ) +

  # Error bars for pooled-era points
  geom_errorbar(
    data = smoker3_pooled_plot,
    aes(
      x = year,
      ymin = estimate - 1.96 * std.error,
      ymax = estimate + 1.96 * std.error
    ),
    width = 0.4,
    linewidth = 0.8,
    color = "#2E8B57"
  ) +

  geom_text(
    data = smoker3_pooled_plot,
    aes(x = year, y = estimate, label = label, vjust = ifelse(year == 2021.5, -1.2, 1.4)),
    color = "#2E8B57",
    size = 3
  ) +

  # Green brackets showing pooled-year ranges
  annotate(
    "segment",
    x = c(2012, 2014, 2020),
    xend = c(2014, 2016, 2023),
    y = y_br,
    yend = y_br,
    linewidth = 1,
    color = "#2E8B57"
  ) +
  annotate(
    "segment",
    x = c(2012, 2014, 2020),
    xend = c(2012, 2014, 2020),
    y = y_br - 0.1,
    yend = y_br + 0.1,
    linewidth = 1,
    color = "#2E8B57"
  ) +
  annotate(
    "segment",
    x = c(2014, 2016, 2023),
    xend = c(2014, 2016, 2023),
    y = y_br - 0.1,
    yend = y_br + 0.1,
    linewidth = 1,
    color = "#2E8B57"
  )

Smoking as a Predictor of Diabetes (BRFSS 2011 - 2023)

Ivan Jimenez- Acosta

2025-12-16