i. Payment penalties

The penalty amounts for each hospital are documented in the CMS FY2022 Final Rule. The file we need is contained within a zip file; we can use this solution for motivation of our approach. First, we can download the file:

# Set location of zip file (publicly available on CMS website)
hrrp_zip <- "https://www.cms.gov/files/zip/fy2022-hrrp-final-rule-supplemental-file.zip"

# Create a temporary file (Credits: rpubs.com/otienodominic/398952)
temp_file <- tempfile()

# Download into the temporary file
download.file(hrrp_zip, temp_file)

Next, we can pluck out the data file we need, save it to our current working directory, and then import it into our session.

# Name of file needed within zip
hrrp_file <- "FY2022_Final_Rule_Supplemental_File.xlsx"

# Unzip, and place the file in the current working directory
unzip(temp_file, hrrp_file, exdir = ".")

# Import the data file into a data frame
hrrp_results <-
  readxl::read_xlsx(
    path = hrrp_file,
    sheet = "FR FY 2022",
    skip = 3,
    na = c("", " ", ".", "-", "NA", "N/A")
  )

# Delete the downloaded file
file.remove(hrrp_file)
unlink(temp_file)

print(hrrp_results, n = 5)

## # A tibble: 3,048 × 36
##   `Hospital CCN` Payme…¹ Payme…² Dual …³ Peer …⁴ Neutr…⁵ Numbe…⁶ ERR f…⁷ Peer …⁸
##   <chr>            <dbl>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl>
## 1 010001           0.990    0.97   0.163       2   0.964     424   1.03    0.996
## 2 010005           0.998    0.15   0.141       2   0.964      36   0.955   0.996
## 3 010006           0.984    1.6    0.107       1   0.964     363   0.961   0.994
## 4 010007           0.995    0.51   0.235       4   0.964       2   0.996   0.998
## 5 010008           1        0      0.227       3   0.964       6   0.998   0.996
## # … with 3,043 more rows, 27 more variables: `Penalty indicator for AMI` <chr>,
## #   `DRG payment ratio for AMI` <dbl>,
## #   `Number of eligible discharges for COPD` <dbl>, `ERR for COPD` <dbl>,
## #   `Peer group median ERR for COPD` <dbl>, `Penalty indicator for COPD` <chr>,
## #   `DRG payment ratio for COPD` <dbl>,
## #   `Number of eligible discharges for HF` <dbl>, `ERR for HF` <dbl>,
## #   `Peer group median ERR for HF` <dbl>, `Penalty indicator for HF` <chr>, …

ii. Hospital information

We need some identifying information for the hospitals such as name and zip code (this is what we’ll use for plotting). This is found on data.cms.gov, which we can import from directly:

# Set path to file (from CMS)
hosp_file <- "https://data.cms.gov/provider-data/sites/default/files/resources/092256becd267d9eeccf73bf7d16c46b_1689206722/Hospital_General_Information.csv"

# Import the file
hosp_info <-
  read_csv(
    file = hosp_file
  )

print(hosp_info, n = 5)

## # A tibble: 5,446 × 38
##   Facili…¹ Facil…² Address City/…³ State ZIP C…⁴ Count…⁵ Telep…⁶ Hospi…⁷ Hospi…⁸
##   <chr>    <chr>   <chr>   <chr>   <chr> <chr>   <chr>   <chr>   <chr>   <chr>  
## 1 010001   SOUTHE… 1108 R… DOTHAN  AL    36301   HOUSTON (334) … Acute … Govern…
## 2 010005   MARSHA… 2505 U… BOAZ    AL    35957   MARSHA… (256) … Acute … Govern…
## 3 010006   NORTH … 1701 V… FLOREN… AL    35630   LAUDER… (256) … Acute … Propri…
## 4 010007   MIZELL… 702 N … OPP     AL    36467   COVING… (334) … Acute … Volunt…
## 5 010008   CRENSH… 101 HO… LUVERNE AL    36049   CRENSH… (334) … Acute … Propri…
## # … with 5,441 more rows, 28 more variables: `Emergency Services` <chr>,
## #   `Meets criteria for promoting interoperability of EHRs` <chr>,
## #   `Hospital overall rating` <chr>, `Hospital overall rating footnote` <dbl>,
## #   `MORT Group Measure Count` <chr>, `Count of Facility MORT Measures` <chr>,
## #   `Count of MORT Measures Better` <chr>,
## #   `Count of MORT Measures No Different` <chr>,
## #   `Count of MORT Measures Worse` <chr>, `MORT Group Footnote` <dbl>, …

iii. Zip code centroids

For simplicity, we are going to use the hospital’s zip code for plotting them on the map. The tigris package provides access to the geometries we need:

# Get set of zip codes for Wisconsin
zips_wi <- tigris::zctas(year = 2010, state = "WI")

A zip code is a region, so we’ll need to reduce those down to a single set of coordinates. To do this, we can use the centroid of the region. The sf package has readily-available functions for this.

# Gather centroids of corrdinates for each zip
zips_wi_centroids <-
  zips_wi %>%
  
  # Get the centroid
  sf::st_centroid() %>%
  
  # Pluck the coordinates
  sf::st_coordinates() %>%
  
  # Make a tibble
  as_tibble() %>%
  
  # Add identifying column
  add_column(
    Zip = zips_wi$ZCTA5CE10
  ) %>%
  
  # Rename columns
  rename(
    lon = X,
    lat = Y
  )

print(zips_wi_centroids, n = 5)

## # A tibble: 774 × 3
##     lon   lat Zip  
##   <dbl> <dbl> <chr>
## 1 -89.8  43.7 53965
## 2 -89.6  44.1 54943
## 3 -89.5  44.2 54966
## 4 -90.7  46.3 54546
## 5 -90.5  46.5 54559
## # … with 769 more rows

iv. Build map data set

Finally, we can clean and combine the previously loaded datasets, extracting the fields we need. First, we’ll filter to Wisconsin hospitals:

map_data <-
  hosp_info %>%
  
  # Filter to Wisconsin hospitals
  filter(
    State == "WI"
  ) %>%
  
  # Keep a few pieces of information
  select(
    FacilityID = `Facility ID`,
    FacilityName = `Facility Name`,
    Zip = `ZIP Code`
  )

Next, we can add the zip code centroid coordinates for the hospital. One additional step we’ll take is to add a small amount of random noise to the coordinates. This will prevent hospitals that share the same zip code from being plotted directly on top of one another.

map_data <-
  map_data %>%
  
  # Join to get the centroid for the hospital's zip code
  inner_join(
    y = zips_wi_centroids,
    by = "Zip"
  ) %>%
  
  # Add random jitter to coordinates
  mutate(
    across(
      c(
        lat,
        lon
      ),
      jitter,
      amount = 0.05
    )
  )

The penalty information can be added next. We’ll add the:

• Penalty amount: The percentage in which reimbursements are reduced by due to excess readmissions
• Dual-eligibility proportion: The percent of discharges in which the patient was dually-eligible for both Medicare and Medicaid benefits
• Peer group: The group that the hospital was compared against to determine penalty status. This is based on the dual-eligibility proportion.
• Cohort-specific penalty indicators: A flag of whether or not the hospital was penalized in each of the six (6) focus cohorts

We’ll do a little clean-up on the cohort indicators for nicer names and presentation.

map_data <-
  map_data %>%
  
  # Join to get HRRP program results
  inner_join(
    y = 
      hrrp_results %>%
      
      # Make cleaner levels for penalty indicators
      mutate(
        across(
          starts_with("Penalty indicator"),
          ~
            case_when(
              .x == "Y" ~ "Yes",
              .x == "N" ~ "No",
              TRUE ~ NA_character_
            )
        )
      ) %>%
      
      # Keep a few columns
      select(
        FacilityID = `Hospital CCN`,
        PeerGroup = `Peer group assignment`,
        DualProportion = `Dual proportion`,
        Penalty = `Payment Reduction Percentage`,
        starts_with("Penalty indicator")
      ) %>%
      
      # Remove the prefix from the cohort columns
      rename_with(
        ~str_remove(.x, "^Penalty indicator for "),
        starts_with("Penalty indicator")
      ),
    by = "FacilityID"
  )

Finally, we’ll add a string that concatenates all of the cohorts that the hospital received penalty for. This will be used as a display for the tooltip in the map.

map_data <-
  map_data %>%
  
  # Join to get list of penalized cohorts (for labeling)
  left_join(
    y = 
      hrrp_results %>%
      
      # Send down the rows
      pivot_longer(
        cols = starts_with("Penalty indicator"),
        names_prefix = "Penalty indicator for "
      ) %>%
      
      # Filter to penalty cohorts
      filter(
        value == "Y"
      ) %>%
      
      # For each hospital
      group_by(`Hospital CCN`) %>%
      
      # Concatenate the list of penalty cohorts
      summarise(
        PenalizedCohortCount = n(),
        PenalizedCohorts = paste(sort(unique(name)), collapse = ", "),
        .groups = "drop"
      ) %>%
      
      # Rename the column
      rename(
        FacilityID = `Hospital CCN`
      ),
    by = "FacilityID"
  ) %>%
  
  # Fill in non-penalized hospitals
  mutate(
    PenalizedCohortCount = coalesce(PenalizedCohortCount, 0),
    PenalizedCohorts = coalesce(PenalizedCohorts, "")
  ) 

print(map_data, n = 5)

## # A tibble: 61 × 16
##   FacilityID Facil…¹ Zip     lon   lat PeerG…² DualP…³ Penalty AMI   COPD  HF   
##   <chr>      <chr>   <chr> <dbl> <dbl>   <dbl>   <dbl>   <dbl> <chr> <chr> <chr>
## 1 520002     ASPIRU… 54481 -89.7  44.6       4   0.262    0.12 No    Yes   Yes  
## 2 520004     MAYO C… 54601 -91.2  43.8       4   0.236    0.81 No    No    No   
## 3 520008     WAUKES… 53188 -88.3  43.0       2   0.139    0.31 Yes   No    Yes  
## 4 520009     ASCENS… 54915 -88.3  44.2       3   0.201    0.08 Yes   No    No   
## 5 520011     MARSHF… 54868 -91.8  45.5       4   0.259    0.37 No    No    Yes  
## # … with 56 more rows, 5 more variables: pneumonia <chr>, CABG <chr>,
## #   `THA/TKA` <chr>, PenalizedCohortCount <dbl>, PenalizedCohorts <chr>, and
## #   abbreviated variable names ¹​FacilityName, ²​PeerGroup, ³​DualProportion