R Notebook

library(tidyverse)
library(survey)
library(ggplot2)
library(plotly)
source("table-builders.R")

acs_data_year = 2024

pus = readRDS(paste0("../../Data/ACS_MEDICAID",acs_data_year,".rds"))

Analysis 1

Source

ACS 2024

Outome

Medicaid, Medical Assistance, or any kind of government-assistance plan for those with low incomes or a disability

Domain of analysis

Children living in the community, with and without disability, who receive medicaid

t1 <- 
  tdev3(data = pus,
       measure = c("POP","MEDICAID"),
       filters = c("POP","AGEP"),
       demographic = "POP == 1 &   
                      AGEP <= 17", 
       strata = "DIS_text",
       strataLevels = c("Yes","No"))

#[[1]]
#[[1]][[1]]

# Population without Disabilities

#[[1]][[1]]

# Population with Disabilities

#[[2]]
#[[2]][[1]]

t1

[[1]]
[[1]][[1]]

[[1]][[2]]

[[2]]
[[2]][[1]]

[[2]][[2]]

NANA

Analysis 2

Source

ACS 2024

Outome

Medicaid, Medical Assistance, or any kind of government-assistance plan for those with low incomes or a disability

Domain of analysis

Adults 18-64 living in the community, with and without disability, who receive medicaid

t2 <- 
  tdev3(data = pus,
       measure = c("POP","MEDICAID"),
       filters = c("POP","AGEP"),
       demographic = "POP == 1 &   
                      AGEP >= 18 &
                      AGEP <= 64", 
       strata = "DIS_text",
       strataLevels = c("Yes","No"))


#[[1]]
#[[1]][[1]]

# Population without Disabilities

#[[1]][[1]]

# Population with Disabilities

#[[2]]
#[[2]][[1]]

t2

[[1]]
[[1]][[1]]

[[1]][[2]]

[[2]]
[[2]][[1]]

[[2]][[2]]

NANA

Medicaid Expansion States

Source: https://www.kff.org/affordable-care-act/state-indicator/state-activity-around-expanding-medicaid-under-the-affordable-care-act/?currentTimeframe=0&sortModel=%7B%22colId%22:%22Location%22,%22sort%22:%22asc%22%7D

expansion = read.csv("../../Data/Medicaid Expansion States/raw_data.csv") %>% 
  filter(!Status.of.Medicaid.Expansion.Decision=="")

expansion <- expansion %>%
  mutate(ExpansionYear = ifelse(is.na(year(as.Date(expansion$Expansion.Implementation.Date, format = "%m/%d/%Y"))),9999,year(as.Date(expansion$Expansion.Implementation.Date, format = "%m/%d/%Y")))) %>%
  select(Location,ExpansionYear)

expansion

NH Yearly Analysis

pnh <- readRDS(paste0("../../Data/ACS_MEDICAID_NH",2008,".rds"))

names(pnh) <- names(pnh) %>% str_replace("pwgtp","PWGTP")

for (i in 2009:2024) {
  print(i)
  tmp <- readRDS(paste0("../../Data/ACS_MEDICAID_NH",i,".rds"))
  names(tmp) <- names(tmp) %>% str_replace("pwgtp","PWGTP")
  pnh <- pnh %>% 
    rbind(tmp)
}

[1] 2009
[1] 2010
[1] 2011
[1] 2012
[1] 2013
[1] 2014
[1] 2015
[1] 2016
[1] 2017
[1] 2018
[1] 2019
[1] 2020
[1] 2021
[1] 2022
[1] 2023
[1] 2024

pops = data.frame()

for (i in 2009:2024) {
  # estimate the pops, prevalence and RR in NH
  out = tdev3(data = pnh %>% filter(YEAR == i),
       measure = c("POP","MEDICAID"),
       filters = c("POP","AGEP"),
       demographic = "POP == 1 &   
                      AGEP >= 18 &
                      AGEP <= 64", 
       strata = "DIS_text",
       strataLevels = c("Yes","No"))
  pops = pops %>% 
    rbind(out[[1]][[1]] %>% mutate(year = i)) %>% rbind(out[[1]][[2]] %>% mutate(year = i))
}

pops <- pops %>% filter(ST_text != "U.S.")

pops

NA
NA

averages <- pops %>%
  mutate(expansionstat = ifelse(year <=2014,"Pre-",
                         ifelse(year > 2014, "Post-",NA)),
         se = (MEDICAID_perc - MEDICAID_perc_LCI)/qnorm(0.975)) %>%
  group_by(`eval(parse(text = strata))`,expansionstat) %>%
  summarise(MEDICAID_perc = mean(MEDICAID_perc),
            se = sqrt(sum(se^2)/n()),
            upperCI= (MEDICAID_perc + qnorm(0.975)*se) %>% round(1),
            lowerCI = (MEDICAID_perc - qnorm(0.975)*se) %>% round(1)) %>%
  filter(expansionstat != "") %>%
  right_join(pops %>%
  mutate(expansionstat = ifelse(year <=2014,"Pre-",
                         ifelse(year > 2014, "Post-",NA))), by = c("eval(parse(text = strata))","expansionstat"))
  

averages %>% 
  arrange(year) 

p <- averages %>% 
  arrange(year) %>%
  ggplot(aes(x =year, group = `eval(parse(text = strata))`)) + 
  geom_ribbon(aes(ymin = lowerCI, ymax = upperCI, fill="steelblue"), alpha = 0.5) + 
  geom_line(aes(y = MEDICAID_perc.x, color = "steelblue", alpha = 0.1)) +
  geom_line(aes(y = MEDICAID_perc.y, color = `eval(parse(text = strata))`)) +
  annotate("text", label="Expansion Implementation Date 8/15/2014", x = 2016.666, y = 20) + 
  geom_vline(xintercept = expansion %>% filter(Location == "New Hampshire") %>% select(ExpansionYear) %>% unlist()) + 
  scale_x_continuous(breaks = 2008:2024) + 
  theme(legend.position = "bottom")

ggplotly(p, width = 1100, height = 1100)%>%
      plotly::layout(legend=list(x=0, 
                                 xanchor='left',
                                 yanchor='bottom',
                                 orientation='h'))