02-analysis-vbac.Rmd

---
title: "AHRQ Vaginal birth after Cesarean"
date: "`r Sys.Date()`"
output:
  html_document:
    df_print: paged
knit: (function(inputFile, encoding) { rmarkdown::render(inputFile, encoding = encoding, output_dir = "docs") })
---

By **Christian McDonald**, Assistant Professor of Practice\
School of Journalism and Media, Moody College of Communication\
University of Texas at Austin

This is an analysis of the Vaginal Birth After Cesareans, as defined by [IQI 22 Vaginal Birth After Cesarean (VBAC) Delivery Rate, Uncomplicated](https://www.qualityindicators.ahrq.gov/Downloads/Modules/IQI/V2020/TechSpecs/IQI_22_Vaginal_Birth_After_Cesarean_(VBAC)_Delivery_Rate_Uncomplicated.pdf) in 2020.

## Definition

### Numerator

Number of vaginal deliveries among cases meeting the inclusion and exclusion rules for the
denominator. Vaginal deliveries are identified by any-listed ICD-10-PCS procedure code for vaginal delivery (VAGDELP* ).

### Denominator

Discharges with an ICD-10-CM diagnosis code for birth delivery outcome (DELOCMD) with any listed ICD-10-CM diagnosis codes for previous Cesarean delivery (PRVBACD).

Exclude cases:

- with any-listed ICD-10-CM diagnosis codes for abnormal presentation, fetal death, or multiple
gestation (Appendix A: PRCSECD )
- with an ungroupable DRG (DRG=999)
- with missing gender (SEX=missing), age (AGE=missing), quarter (DQTR=missing), year
(YEAR=missing) or principal diagnosis (DX1=missing)

```{r setup, echo=T, results='hide', message=F, warning=F}
library(tidyverse)
library(janitor)
library(DT)
library(tigris)

# suppresses grouping warning
options(dplyr.summarise.inform = FALSE)
```

## Import of deliveries

I start here with "All deliveries, identified by any-listed ICD-10-CM diagnosis code for outcome of delivery (DELOCMD)". This was processed in `01-process-loop`.

```{r imports}
# set test flag to FALSE to run production data
test_flag <- F

### test data
path_test <- "data-test/ahrq_del_all_loop_test.rds"

### production data
path_prod <- "data-processed/ahrq_del_cleaned.rds"

### import based on flag
if (test_flag == T) del <- read_rds(path_test) else del <- read_rds(path_prod)

del %>% nrow()
```

## Set up various processing lists

These are lists from the IQI 22 reference used for filtering and such. See 01-process-lists for definitions.

```{r lists}
# diagnostic columns
diag_cols <- read_rds("procedures-lists/cols_diag.rds") %>% .$diag

# surgical procedure columns
surg_cols <- read_rds("procedures-lists/cols_surg.rds") %>% .$surg

# appendix a complications list
prcsecd_list <- read_rds("procedures-lists/ahrq_prcsecd.rds") %>% .$prcsecd

# Previous Cesarean delivery diagnosis codes: (PRVBACD)
prvbacd_list <- read_rds("procedures-lists/ahrq_prvbacd.rds") %>% .$prvbacd

# Vaginal deliveries (VAGDELP)
vagdelp_list <- read_rds("procedures-lists/ahrq_vagdelp.rds") %>% .$vagdelp
```

## Filter deliveries to set denominator

### Filter out complicated deliveries

Start by filtering out rows "with any-listed ICD-10-CM diagnosis codes for abnormal presentation, preterm, fetal death, or multiple gestation (Appendix A: PRCSECD)".


```{r prcsecd}
del_ucmp <- del %>% 
  filter_at(
    vars(all_of(diag_cols)),
    all_vars(
      !(. %in% prcsecd_list)
    )
  )

del_ucmp %>% nrow()
```

### Filter out ungroupable DRG (DRG=999)

There are two Diagnosis Related Groups in the data, `MS_DRG`: Centers for Medicare and Medicaid Services (CMS) Diagnosis Related Group (DRG), as assigned for hospital payment for Medicare beneficiaries; and `APR_DRG`: All Patient Refined (APR) Diagnosis Related Group (DRG) as assigned by 3M APR-DRG Grouper). I filter out records with "999" in either of them.

```{r drg_ungroupable}
del_ucmp <- del_ucmp %>% 
  filter(
    (APR_DRG != "999"),
    (MS_DRG != "999")
  )

del_ucmp %>% nrow()
```

### Filter to keep only previous Cesarean delivery (PRVBACD).

```{r prvbacd}
del_denom <- del_ucmp %>% 
  filter_at(
    vars(all_of(diag_cols)),
    any_vars(
      . %in% prvbacd_list
    )
  )

del_denom %>% nrow()
```

This completes our "denominator" set of data. This is saved in `del_ucmp`.

## Numerator updates

Number of vaginal deliveries among cases meeting the inclusion and exclusion rules for the denominator. Vaginal deliveries are identified by any-listed ICD-10-PCS procedure code for vaginal delivery (VAGDELP).

We'll add a new column `VAGDELP` based on the codes. If the code is present, the value is `TRUE`.

```{r vagdelp}
ahrq_vbac <- del_denom %>% 
  mutate(
    VBAC = case_when(
      PRINC_SURG_PROC_CODE %in% vagdelp_list ~ TRUE,
      OTH_SURG_PROC_CODE_1 %in% vagdelp_list ~ TRUE,
      OTH_SURG_PROC_CODE_2 %in% vagdelp_list ~ TRUE,
      OTH_SURG_PROC_CODE_3 %in% vagdelp_list ~ TRUE,
      OTH_SURG_PROC_CODE_4 %in% vagdelp_list ~ TRUE,
      OTH_SURG_PROC_CODE_5 %in% vagdelp_list ~ TRUE,
      OTH_SURG_PROC_CODE_6 %in% vagdelp_list ~ TRUE,
      OTH_SURG_PROC_CODE_7 %in% vagdelp_list ~ TRUE,
      OTH_SURG_PROC_CODE_8 %in% vagdelp_list ~ TRUE,
      OTH_SURG_PROC_CODE_9 %in% vagdelp_list ~ TRUE,
      OTH_SURG_PROC_CODE_10 %in% vagdelp_list ~ TRUE,
      OTH_SURG_PROC_CODE_11 %in% vagdelp_list ~ TRUE,
      OTH_SURG_PROC_CODE_12 %in% vagdelp_list ~ TRUE,
      OTH_SURG_PROC_CODE_13 %in% vagdelp_list ~ TRUE,
      OTH_SURG_PROC_CODE_14 %in% vagdelp_list ~ TRUE,
      OTH_SURG_PROC_CODE_15 %in% vagdelp_list ~ TRUE,
      OTH_SURG_PROC_CODE_16 %in% vagdelp_list ~ TRUE,
      OTH_SURG_PROC_CODE_17 %in% vagdelp_list ~ TRUE,
      OTH_SURG_PROC_CODE_18 %in% vagdelp_list ~ TRUE,
      OTH_SURG_PROC_CODE_19 %in% vagdelp_list ~ TRUE,
      OTH_SURG_PROC_CODE_20 %in% vagdelp_list ~ TRUE,
      OTH_SURG_PROC_CODE_21 %in% vagdelp_list ~ TRUE,
      OTH_SURG_PROC_CODE_22 %in% vagdelp_list ~ TRUE,
      OTH_SURG_PROC_CODE_23 %in% vagdelp_list ~ TRUE,
      OTH_SURG_PROC_CODE_24 %in% vagdelp_list ~ TRUE,
      TRUE                     ~  FALSE
    )
  )

ahrq_vbac %>% nrow()
```

This completes the processing of the data. `ahrq_vbac` is our working dataframe.

---

## VBAC analysis

- Do hospitals perform them?
- How many performed in 2019?

### VBAC by hospital by year

This data is filtered for hospitals that have 30+ deliveries of patients with previous ceseareans in a given year.

```{r ahrq_vbac_rate_hosp_yr}
ahrq_vbac_rate_hosp_yr <- ahrq_vbac %>% 
  group_by(YR, THCIC_ID, PROVIDER_NAME) %>% 
  count(VBAC) %>% 
  # pivot to get rate for vbac
  pivot_wider(names_from = VBAC, values_from = n) %>% 
  rename(
    NVBAC = `FALSE`,
    YVBAC = `TRUE`
  ) %>% 
  mutate(
    PCTOT = NVBAC + YVBAC,
    VBACRATE = ((YVBAC / PCTOT) * 100) %>% round_half_up(1)
  ) %>% 
  filter(PCTOT >= 30 ) %>% 
  arrange(YR, PROVIDER_NAME)
  
ahrq_vbac_rate_hosp_yr %>% datatable()
```

### Rates in a table by year

```{r vbac_hosp_yr_table}
vbac_hosp_yr_table <- ahrq_vbac_rate_hosp_yr %>%
  select(YR, THCIC_ID, PROVIDER_NAME, VBACRATE) %>% 
  pivot_wider(names_from = YR, values_from = c(VBACRATE)) %>% 
  arrange(PROVIDER_NAME)

vbac_hosp_yr_table %>% datatable()
```

## VBACs by Laredo hospitals

```{r}
vbac_hosp_yr_table %>% 
  filter(str_detect(PROVIDER_NAME, "Laredo"))
```
## Summaries

### Statewide VBAC percentage across data

```{r ahrq_vbac_rate_tx}
ahrq_vbac %>% 
  tabyl(VBAC) %>% 
  rename(count = n) %>% 
  adorn_pct_formatting()
```

```{r ahrq_vbac_rate_tx_summary}
ahrq_vbac_rate_tx_summary <- ahrq_vbac %>% 
  group_by(VBAC) %>% 
  summarize(CNT = n()) %>% 
  pivot_wider(names_from = VBAC, values_from = CNT) %>% 
  rename(
    NVBAC_CNT = "FALSE",
    VBAC_CNT = "TRUE"
  ) %>% 
  mutate(
    TOTAL = NVBAC_CNT + VBAC_CNT
  ) %>% 
  mutate(
   SUMMARY = "TX",
   CATEGORY = "VAGINAL_BIRTH_AFTER_CESAREAN",
   MEASUREMENT = "RATE",
   VALUE = round_half_up((VBAC_CNT / TOTAL) * 100,1) # RATE
  ) %>%
  select(SUMMARY, CATEGORY, MEASUREMENT, VALUE)

ahrq_vbac_rate_tx_summary
```

### VBAC rate statewide by year

Excludes hospitals with fewer than 30 deliveries a year.

```{r ahrq_vbac_rate_tx_yr_summary}
ahrq_vbac_rate_tx_yr_summary <- ahrq_vbac %>% 
  group_by(VBAC, YR) %>% 
  summarize(CNT = n()) %>% 
  pivot_wider(names_from = VBAC, values_from = CNT) %>% 
  rename(
    NVBAC_CNT = "FALSE",
    VBAC_CNT = "TRUE"
  ) %>% 
  mutate(
    TOTAL = NVBAC_CNT + VBAC_CNT
  ) %>% 
  filter(TOTAL >= 30) %>% 
  mutate(
   SUMMARY = "TX",
   CATEGORY = "VAGINAL_BIRTH_AFTER_CESAREAN",
   MEASUREMENT = "RATE",
   VALUE = round_half_up((VBAC_CNT / TOTAL) * 100,1) # RATE
  ) %>%
  select(YR, SUMMARY, CATEGORY, MEASUREMENT, VALUE)

ahrq_vbac_rate_tx_yr_summary
```

### VBAC rate by hospital: Averaged, by year

```{r ahrq_vbac_rate_hosp_yr_summary}
ahrq_vbac_rate_hosp_yr_summary <- ahrq_vbac_rate_hosp_yr %>%
  ungroup() %>% 
  group_by(YR) %>% 
  summarize(
    SUMMARY = "HOSPITAL",
    CATEGORY = "VAGINAL_BIRTH_AFTER_CESAREAN",
    MEASUREMENT = "MEAN_OF_RATE",
    # MEDIAN = median(PCRATE),
    VALUE = round_half_up(mean(VBACRATE),1) #MEAN
  )

ahrq_vbac_rate_hosp_yr_summary
```

## Exports

Writing out some summary data for other notebooks, including:

- ahrq_vbac_rate_hosp_yr
- ahrq_vbac_rate_hosp_yr_summary
- ahrq_vbac_rate_tx_yr_summary

```{r exports}
if (test_flag == F) ahrq_vbac_rate_tx_yr_summary %>%
  write_rds("data-processed/ahrq_vbac_rate_tx_yr_summary.rds")

if (test_flag == F) ahrq_vbac_rate_hosp_yr_summary %>%
  write_rds("data-processed/ahrq_vbac_rate_hosp_yr_summary.rds")

if (test_flag == F) ahrq_vbac_rate_hosp_yr %>%
  write_rds("data-processed/ahrq_vbac_rate_hosp_yr.rds")
```


## Closing

```{r end}
beepr::beep(4)
```