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 — CHECKUP1
# -------------------------
tx_checkup1_2011 <- BRFSS2011 %>%
filter(x.state == 48) %>%
select(x.psu, x.ststr, x.llcpwt, diabete3, checkup1) %>%
mutate(
diabete3 = ifelse(diabete3 %in% c(7, 9), NA, diabete3),
diabetes_dx = case_when(diabete3 == 1 ~ 1, diabete3 == 3 ~ 0, TRUE ~ NA_real_),
checkup1 = ifelse(checkup1 %in% c(7, 9), NA, checkup1)
) %>%
filter(!is.na(diabetes_dx), !is.na(checkup1)) %>%
mutate(checkup1 = relevel(factor(checkup1, levels = c(1,2,3,4,8)), ref = "1"))
des_checkup1_2011 <- svydesign(ids=~x.psu, strata=~x.ststr, weights=~x.llcpwt,
data=tx_checkup1_2011, nest=TRUE)
out_checkup1_2011 <- broom::tidy(
svyglm(diabetes_dx ~ checkup1, design=des_checkup1_2011, family=quasibinomial()),
conf.int=FALSE
) %>%
filter(term != "(Intercept)") %>%
mutate(year=2011, predictor="checkup1")
# -------------------------
# 2012 — CHECKUP1
# -------------------------
tx_checkup1_2012 <- BRFSS2012 %>%
filter(x.state == 48) %>%
select(x.psu, x.ststr, x.llcpwt, diabete3, checkup1) %>%
mutate(
diabete3 = ifelse(diabete3 %in% c(7, 9), NA, diabete3),
diabetes_dx = case_when(diabete3 == 1 ~ 1, diabete3 == 3 ~ 0, TRUE ~ NA_real_),
checkup1 = ifelse(checkup1 %in% c(7, 9), NA, checkup1)
) %>%
filter(!is.na(diabetes_dx), !is.na(checkup1)) %>%
mutate(checkup1 = relevel(factor(checkup1, levels = c(1,2,3,4,8)), ref = "1"))
des_checkup1_2012 <- svydesign(ids=~x.psu, strata=~x.ststr, weights=~x.llcpwt,
data=tx_checkup1_2012, nest=TRUE)
out_checkup1_2012 <- broom::tidy(
svyglm(diabetes_dx ~ checkup1, design=des_checkup1_2012, family=quasibinomial()),
conf.int=FALSE
) %>%
filter(term != "(Intercept)") %>%
mutate(year=2012, predictor="checkup1")
# -------------------------
# 2013 — CHECKUP1
# -------------------------
tx_checkup1_2013 <- BRFSS2013 %>%
filter(x.state == 48) %>%
select(x.psu, x.ststr, x.llcpwt, diabete3, checkup1) %>%
mutate(
diabete3 = ifelse(diabete3 %in% c(7, 9), NA, diabete3),
diabetes_dx = case_when(diabete3 == 1 ~ 1, diabete3 == 3 ~ 0, TRUE ~ NA_real_),
checkup1 = ifelse(checkup1 %in% c(7, 9), NA, checkup1)
) %>%
filter(!is.na(diabetes_dx), !is.na(checkup1)) %>%
mutate(checkup1 = relevel(factor(checkup1, levels = c(1,2,3,4,8)), ref = "1"))
des_checkup1_2013 <- svydesign(ids=~x.psu, strata=~x.ststr, weights=~x.llcpwt,
data=tx_checkup1_2013, nest=TRUE)
out_checkup1_2013 <- broom::tidy(
svyglm(diabetes_dx ~ checkup1, design=des_checkup1_2013, family=quasibinomial()),
conf.int=FALSE
) %>%
filter(term != "(Intercept)") %>%
mutate(year=2013, predictor="checkup1")
# -------------------------
# 2014 — CHECKUP1
# -------------------------
tx_checkup1_2014 <- BRFSS2014 %>%
filter(x.state == 48) %>%
select(x.psu, x.ststr, x.llcpwt, diabete3, checkup1) %>%
mutate(
diabete3 = ifelse(diabete3 %in% c(7, 9), NA, diabete3),
diabetes_dx = case_when(diabete3 == 1 ~ 1, diabete3 == 3 ~ 0, TRUE ~ NA_real_),
checkup1 = ifelse(checkup1 %in% c(7, 9), NA, checkup1)
) %>%
filter(!is.na(diabetes_dx), !is.na(checkup1)) %>%
mutate(checkup1 = relevel(factor(checkup1, levels = c(1,2,3,4,8)), ref = "1"))
des_checkup1_2014 <- svydesign(ids=~x.psu, strata=~x.ststr, weights=~x.llcpwt,
data=tx_checkup1_2014, nest=TRUE)
out_checkup1_2014 <- broom::tidy(
svyglm(diabetes_dx ~ checkup1, design=des_checkup1_2014, family=quasibinomial()),
conf.int=FALSE
) %>%
filter(term != "(Intercept)") %>%
mutate(year=2014, predictor="checkup1")
# -------------------------
# 2015 — CHECKUP1
# -------------------------
tx_checkup1_2015 <- BRFSS2015 %>%
filter(x.state == 48) %>%
select(x.psu, x.ststr, x.llcpwt, diabete3, checkup1) %>%
mutate(
diabete3 = ifelse(diabete3 %in% c(7, 9), NA, diabete3),
diabetes_dx = case_when(diabete3 == 1 ~ 1, diabete3 == 3 ~ 0, TRUE ~ NA_real_),
checkup1 = ifelse(checkup1 %in% c(7, 9), NA, checkup1)
) %>%
filter(!is.na(diabetes_dx), !is.na(checkup1)) %>%
mutate(checkup1 = relevel(factor(checkup1, levels = c(1,2,3,4,8)), ref = "1"))
des_checkup1_2015 <- svydesign(ids=~x.psu, strata=~x.ststr, weights=~x.llcpwt,
data=tx_checkup1_2015, nest=TRUE)
out_checkup1_2015 <- broom::tidy(
svyglm(diabetes_dx ~ checkup1, design=des_checkup1_2015, family=quasibinomial()),
conf.int=FALSE
) %>%
filter(term != "(Intercept)") %>%
mutate(year=2015, predictor="checkup1")
# -------------------------
# 2016 — CHECKUP1
# -------------------------
tx_checkup1_2016 <- BRFSS2016 %>%
filter(x.state == 48) %>%
select(x.psu, x.ststr, x.llcpwt, diabete3, checkup1) %>%
mutate(
diabete3 = ifelse(diabete3 %in% c(7, 9), NA, diabete3),
diabetes_dx = case_when(diabete3 == 1 ~ 1, diabete3 == 3 ~ 0, TRUE ~ NA_real_),
checkup1 = ifelse(checkup1 %in% c(7, 9), NA, checkup1)
) %>%
filter(!is.na(diabetes_dx), !is.na(checkup1)) %>%
mutate(checkup1 = relevel(factor(checkup1, levels = c(1,2,3,4,8)), ref = "1"))
des_checkup1_2016 <- svydesign(ids=~x.psu, strata=~x.ststr, weights=~x.llcpwt,
data=tx_checkup1_2016, nest=TRUE)
out_checkup1_2016 <- broom::tidy(
svyglm(diabetes_dx ~ checkup1, design=des_checkup1_2016, family=quasibinomial()),
conf.int=FALSE
) %>%
filter(term != "(Intercept)") %>%
mutate(year=2016, predictor="checkup1")
# -------------------------
# 2017 — CHECKUP1
# -------------------------
tx_checkup1_2017 <- BRFSS2017 %>%
filter(x.state == 48) %>%
select(x.psu, x.ststr, x.llcpwt, diabete3, checkup1) %>%
mutate(
diabete3 = ifelse(diabete3 %in% c(7, 9), NA, diabete3),
diabetes_dx = case_when(diabete3 == 1 ~ 1, diabete3 == 3 ~ 0, TRUE ~ NA_real_),
checkup1 = ifelse(checkup1 %in% c(7, 9), NA, checkup1)
) %>%
filter(!is.na(diabetes_dx), !is.na(checkup1)) %>%
mutate(checkup1 = relevel(factor(checkup1, levels = c(1,2,3,4,8)), ref = "1"))
des_checkup1_2017 <- svydesign(ids=~x.psu, strata=~x.ststr, weights=~x.llcpwt,
data=tx_checkup1_2017, nest=TRUE)
out_checkup1_2017 <- broom::tidy(
svyglm(diabetes_dx ~ checkup1, design=des_checkup1_2017, family=quasibinomial()),
conf.int=FALSE
) %>%
filter(term != "(Intercept)") %>%
mutate(year=2017, predictor="checkup1")
# -------------------------
# 2018 — CHECKUP1
# -------------------------
tx_checkup1_2018 <- BRFSS2018 %>%
filter(x.state == 48) %>%
select(x.psu, x.ststr, x.llcpwt, diabete3, checkup1) %>%
mutate(
diabete3 = ifelse(diabete3 %in% c(7, 9), NA, diabete3),
diabetes_dx = case_when(diabete3 == 1 ~ 1, diabete3 == 3 ~ 0, TRUE ~ NA_real_),
checkup1 = ifelse(checkup1 %in% c(7, 9), NA, checkup1)
) %>%
filter(!is.na(diabetes_dx), !is.na(checkup1)) %>%
mutate(checkup1 = relevel(factor(checkup1, levels = c(1,2,3,4,8)), ref = "1"))
des_checkup1_2018 <- svydesign(ids=~x.psu, strata=~x.ststr, weights=~x.llcpwt,
data=tx_checkup1_2018, nest=TRUE)
out_checkup1_2018 <- broom::tidy(
svyglm(diabetes_dx ~ checkup1, design=des_checkup1_2018, family=quasibinomial()),
conf.int=FALSE
) %>%
filter(term != "(Intercept)") %>%
mutate(year=2018, predictor="checkup1")
# -------------------------
# 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 — CHECKUP1
# -------------------------
tx_checkup1_2019 <- BRFSS2019 %>%
filter(x.state == 48) %>%
select(x.psu, x.ststr, x.llcpwt, diabete4, checkup1) %>%
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_),
checkup1 = ifelse(checkup1 %in% c(7, 9), NA, checkup1)
) %>%
filter(!is.na(diabetes_dx), !is.na(checkup1)) %>%
mutate(checkup1 = relevel(factor(checkup1, levels = c(1,2,3,4,8)), ref = "1"))
des_checkup1_2019 <- svydesign(ids=~x.psu, strata=~x.ststr, weights=~x.llcpwt,
data=tx_checkup1_2019, nest=TRUE)
out_checkup1_2019 <- broom::tidy(
svyglm(diabetes_dx ~ checkup1, design=des_checkup1_2019, family=quasibinomial()),
conf.int=FALSE
) %>%
filter(term != "(Intercept)") %>%
mutate(year=2019, predictor="checkup1")
# -------------------------
# NOTE:
# BRFSS 2020 Texas data contain five singleton strata, each with exactly
# one PSU and one observation. This violates svyglm() requirements and
# produces the error: "Stratum (...) has only one PSU at stage 1".
# To proceed, these five observations (one per stratum) are dropped from
# the 2020 analytic dataset ONLY.
# Dropped strata: 481012, 481052, 481062, 481092, 481172
# 2020 — CHECKUP1
# -------------------------
tx_checkup1_2020 <- BRFSS2020 %>%
filter(x.state == 48) %>%
select(x.psu, x.ststr, x.llcpwt, diabete4, checkup1) %>%
rename(diabete3 = diabete4) %>% # keep diabetes name consistent in code
mutate(
diabete3 = ifelse(diabete3 %in% c(7, 9), NA, diabete3),
diabetes_dx = case_when(diabete3 == 1 ~ 1, diabete3 == 3 ~ 0, TRUE ~ NA_real_),
checkup1 = ifelse(checkup1 %in% c(7, 9), NA, checkup1)
) %>%
filter(!is.na(diabetes_dx), !is.na(checkup1)) %>%
mutate(
checkup1 = relevel(factor(checkup1, levels = c(1,2,3,4,8)), ref = "1") # ref = within past year
)
# Drop singleton strata that trigger "only one PSU" error
drop_strata_2020 <- c(481012, 481052, 481062, 481092, 481172)
tx_checkup1_2020 <- tx_checkup1_2020 %>%
filter(!(as.numeric(x.ststr) %in% drop_strata_2020))
des_checkup1_2020 <- svydesign(
ids = ~x.psu,
strata = ~x.ststr,
weights = ~x.llcpwt,
data = tx_checkup1_2020,
nest = TRUE
)
out_checkup1_2020 <- broom::tidy(
svyglm(diabetes_dx ~ checkup1, design = des_checkup1_2020, family = quasibinomial()),
conf.int = FALSE
) %>%
filter(term != "(Intercept)") %>%
mutate(year = 2020, predictor = "checkup1")
# -------------------------
# NOTE:
# BRFSS 2021 Texas data contain four singleton strata, each with exactly
# one PSU and one observation. This violates svyglm() requirements and
# produces the error: "Stratum (...) has only one PSU at stage 1".
# To proceed, these four observations (one per stratum) are dropped from
# the 2021 analytic dataset ONLY.
# Dropped strata: 482099, 482149, 482179, 482189
# 2021 — CHECKUP1
# -------------------------
tx_checkup1_2021 <- BRFSS2021 %>%
filter(x.state == 48) %>%
select(x.psu, x.ststr, x.llcpwt, diabete4, checkup1) %>%
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_),
checkup1 = ifelse(checkup1 %in% c(7, 9), NA, checkup1)
) %>%
filter(!is.na(diabetes_dx), !is.na(checkup1)) %>%
mutate(
checkup1 = relevel(factor(checkup1, levels = c(1,2,3,4,7,8,9)), ref = "1")
) %>%
filter(!x.ststr %in% c(482099, 482149, 482179, 482189))
des_checkup1_2021 <- svydesign(
ids = ~x.psu,
strata = ~x.ststr,
weights = ~x.llcpwt,
data = tx_checkup1_2021,
nest = TRUE
)
out_checkup1_2021 <- broom::tidy(
svyglm(diabetes_dx ~ checkup1, design = des_checkup1_2021, family = quasibinomial()),
conf.int = FALSE
) %>%
filter(term != "(Intercept)") %>%
mutate(year = 2021, predictor = "checkup1")
# -------------------------
# NOTE:
# BRFSS 2022 Texas data contain eight singleton strata, each with exactly
# one PSU and one observation. This violates svyglm() requirements and
# produces the error: "Stratum (...) has only one PSU at stage 1".
# To proceed, these eight observations (one per stratum) are dropped from
# the 2022 analytic dataset ONLY.
# Dropped strata: 482029, 482039, 482119, 482139, 482199, 482219, 482229, 482259
# 2022 — CHECKUP1
# -------------------------
tx_checkup1_2022 <- BRFSS2022 %>%
filter(x.state == 48) %>%
select(x.psu, x.ststr, x.llcpwt, diabete4, checkup1) %>%
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_),
checkup1 = ifelse(checkup1 %in% c(7, 9), NA, checkup1)
) %>%
filter(!is.na(diabetes_dx), !is.na(checkup1)) %>%
mutate(
checkup1 = relevel(factor(checkup1, levels = c(1,2,3,4,7,8,9)), ref = "1")
) %>%
filter(!x.ststr %in% c(482029, 482039, 482119, 482139, 482199, 482219, 482229, 482259))
des_checkup1_2022 <- svydesign(
ids = ~x.psu,
strata = ~x.ststr,
weights = ~x.llcpwt,
data = tx_checkup1_2022,
nest = TRUE
)
out_checkup1_2022 <- broom::tidy(
svyglm(diabetes_dx ~ checkup1, design = des_checkup1_2022, family = quasibinomial()),
conf.int = FALSE
) %>%
filter(term != "(Intercept)") %>%
mutate(year = 2022, predictor = "checkup1")
# -------------------------
# NOTE:
# BRFSS 2023 TX contains three strata with only one PSU at stage 1, which causes
# svydesign/svyglm to error: "Stratum (...) has only one PSU at stage 1".
# Two of these strata (482149, 482169) are singleton in the raw TX data,
# and one stratum (481102) becomes singleton after analytic filtering.
# Each stratum contains exactly one observation.
# To proceed, these 3 observations are dropped from the 2023 analytic dataset ONLY:
# 481102, 482149, 482169.
# 2023 — CHECKUP1
# -------------------------
tx_checkup1_2023 <- BRFSS2023 %>%
filter(x.state == 48) %>%
select(x.psu, x.ststr, x.llcpwt, diabete4, checkup1) %>%
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_),
checkup1 = ifelse(checkup1 %in% c(7, 9), NA, checkup1)
) %>%
filter(!is.na(diabetes_dx), !is.na(checkup1)) %>%
mutate(
checkup1 = relevel(factor(checkup1, levels = c(1,2,3,4,8)), ref = "1")
) %>%
# FIX: Drop singleton-PSU strata that cause svydesign/svyglm to error in 2023
filter(!x.ststr %in% c(481102, 482149, 482169))
des_checkup1_2023 <- svydesign(
ids = ~x.psu,
strata = ~x.ststr,
weights = ~x.llcpwt,
data = tx_checkup1_2023,
nest = TRUE
)
out_checkup1_2023 <- broom::tidy(
svyglm(diabetes_dx ~ checkup1, design = des_checkup1_2023, family = quasibinomial()),
conf.int = FALSE
) %>%
filter(term != "(Intercept)") %>%
mutate(year = 2023, predictor = "checkup1")
# Pool 2012–2014 (pre-ACA)
pool_checkup1_2012_2014 <- bind_rows(
tx_checkup1_2012 %>% select(x.psu,x.ststr,x.llcpwt,diabetes_dx,checkup1) %>%
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), checkup1=as.character(checkup1)
),
tx_checkup1_2013 %>% select(x.psu,x.ststr,x.llcpwt,diabetes_dx,checkup1) %>%
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), checkup1=as.character(checkup1)
),
tx_checkup1_2014 %>% select(x.psu,x.ststr,x.llcpwt,diabetes_dx,checkup1) %>%
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), checkup1=as.character(checkup1)
)
) %>%
mutate(
checkup1 = relevel(factor(checkup1, levels = c("1","2","3","4","8")), ref = "1")
) # baseline is 1 by default
des_checkup1_2012_2014 <- svydesign(ids=~x.psu, strata=~x.ststr, weights=~x.llcpwt,
data=pool_checkup1_2012_2014, nest=TRUE)
out_checkup1_2012_2014 <- broom::tidy(
svyglm(diabetes_dx ~ checkup1, design=des_checkup1_2012_2014, family=quasibinomial()),
conf.int=FALSE
)
strong_checkup1_2012_2014 <- out_checkup1_2012_2014 %>%
filter(term!="(Intercept)") %>%
slice_min(p.value, n=1, with_ties=FALSE)
strong_checkup1_2012_2014
# Pool 2014–2016 (early ACA)
pool_checkup1_2014_2016 <- bind_rows(
tx_checkup1_2014 %>% select(x.psu,x.ststr,x.llcpwt,diabetes_dx,checkup1) %>%
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), checkup1=as.character(checkup1)
),
tx_checkup1_2015 %>% select(x.psu,x.ststr,x.llcpwt,diabetes_dx,checkup1) %>%
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), checkup1=as.character(checkup1)
),
tx_checkup1_2016 %>% select(x.psu,x.ststr,x.llcpwt,diabetes_dx,checkup1) %>%
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), checkup1=as.character(checkup1)
)
) %>%
mutate(
checkup1 = relevel(factor(checkup1, levels = c("1","2","3","4","8")), ref = "1")
) # baseline is 1 by default
des_checkup1_2014_2016 <- svydesign(ids=~x.psu, strata=~x.ststr, weights=~x.llcpwt,
data=pool_checkup1_2014_2016, nest=TRUE)
out_checkup1_2014_2016 <- broom::tidy(
svyglm(diabetes_dx ~ checkup1, design=des_checkup1_2014_2016, family=quasibinomial()),
conf.int=FALSE
)
strong_checkup1_2014_2016 <- out_checkup1_2014_2016 %>%
filter(term!="(Intercept)") %>%
slice_min(p.value, n=1, with_ties=FALSE)
strong_checkup1_2014_2016
# Pool 2020–2023 (COVID era)
pool_checkup1_2020_2023 <- bind_rows(
tx_checkup1_2020 %>% select(x.psu,x.ststr,x.llcpwt,diabetes_dx,checkup1) %>%
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), checkup1=as.character(checkup1)
),
tx_checkup1_2021 %>% select(x.psu,x.ststr,x.llcpwt,diabetes_dx,checkup1) %>%
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), checkup1=as.character(checkup1)
),
tx_checkup1_2022 %>% select(x.psu,x.ststr,x.llcpwt,diabetes_dx,checkup1) %>%
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), checkup1=as.character(checkup1)
),
tx_checkup1_2023 %>% select(x.psu,x.ststr,x.llcpwt,diabetes_dx,checkup1) %>%
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), checkup1=as.character(checkup1)
)
) %>%
mutate(
checkup1 = relevel(factor(checkup1, levels = c("1","2","3","4","8")), ref = "1")
) # baseline is 1 by default
des_checkup1_2020_2023 <- svydesign(ids=~x.psu, strata=~x.ststr, weights=~x.llcpwt,
data=pool_checkup1_2020_2023, nest=TRUE)
out_checkup1_2020_2023 <- broom::tidy(
svyglm(diabetes_dx ~ checkup1, design=des_checkup1_2020_2023, family=quasibinomial()),
conf.int=FALSE
)
strong_checkup1_2020_2023 <- out_checkup1_2020_2023 %>%
filter(term!="(Intercept)") %>%
slice_min(p.value, n=1, with_ties=FALSE)
strong_checkup1_2020_2023
# Combine all years
checkup1_all <- bind_rows(
out_checkup1_2011 %>% mutate(year = 2011),
out_checkup1_2012 %>% mutate(year = 2012),
out_checkup1_2013 %>% mutate(year = 2013),
out_checkup1_2014 %>% mutate(year = 2014),
out_checkup1_2015 %>% mutate(year = 2015),
out_checkup1_2016 %>% mutate(year = 2016),
out_checkup1_2017 %>% mutate(year = 2017),
out_checkup1_2018 %>% mutate(year = 2018),
out_checkup1_2019 %>% mutate(year = 2019),
out_checkup1_2020 %>% mutate(year = 2020),
out_checkup1_2021 %>% mutate(year = 2021),
out_checkup1_2022 %>% mutate(year = 2022),
out_checkup1_2023 %>% mutate(year = 2023)
)
# Keep only the CHECKUP1 coefficients (drop intercept)
checkup1_strongest <- checkup1_all %>%
filter(term != "(Intercept)") %>%
group_by(year) %>%
slice_min(order_by = p.value, n = 1, with_ties = FALSE) %>%
ungroup()
# Build checkup1_period_table from pooled strongest objects
checkup1_period_table <- bind_rows(
strong_checkup1_2012_2014 %>% mutate(period = "2012–2014"),
strong_checkup1_2014_2016 %>% mutate(period = "2014–2016"),
strong_checkup1_2020_2023 %>% mutate(period = "2020–2023")
)
# Pooled-period points (from checkup1_period_table)
checkup1_pooled_plot <- checkup1_period_table %>%
mutate(
year = c(2013, 2015, 2021.5),
label = term
)
# Brackets below lowest CI
y_ci_min <- min(checkup1_strongest$estimate - 1.96 * checkup1_strongest$std.error, na.rm = TRUE)
y_ci_max <- max(checkup1_strongest$estimate + 1.96 * checkup1_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(checkup1_strongest, aes(x = year, y = estimate)) +
geom_line() +
geom_point(size = 2.5) +
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 = 3,
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(checkup1_strongest$estimate, na.rm = TRUE),
label = "ACA",
color = "#1f77b4",
angle = 90,
hjust = 5,
vjust = -0.6,
size = 4
) +
annotate(
"text",
x = 2020,
y = max(checkup1_strongest$estimate, na.rm = TRUE),
label = "COVID",
color = "#d62728",
angle = 90,
hjust = 3.5,
vjust = -0.6,
size = 4
) +
scale_x_continuous(breaks = 2011:2023) +
labs(
x = "Year",
y = "Log-odds coefficient",
title = "Most Statistically Significant Checkup 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")
) +
geom_point(
data = checkup1_pooled_plot,
aes(x = year, y = estimate),
color = "#2E8B57",
size = 2.5
) +
geom_errorbar(
data = checkup1_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 = checkup1_pooled_plot,
aes(x = year, y = estimate, label = label,
vjust = ifelse(year == 2021.5, -1.2, 1.4)),
color = "#2E8B57",
size = 3
) +
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.2,
yend = y_br + 0.2,
linewidth = 1,
color = "#2E8B57"
) +
annotate(
"segment",
x = c(2014, 2016, 2023),
xend = c(2014, 2016, 2023),
y = y_br - 0.2,
yend = y_br + 0.2,
linewidth = 1,
color = "#2E8B57"
)