Building a filterable map with {leaflet} and {crosstalk} for hospital readmissions

Healthcare

Maps

Javascript

For runtime-free analytical content delivery

Author

Alex Zajichek

Published

June 18, 2022

Want to make this map widget? Read On! 👇👇👇

Hospital Readmissions Reduction Program (HRRP) FY 2022 - WI Map

Payment Reduction

Cohort-specific penalty:

AMI

CABG

COPD

PNEU

THA/TKA

So you run R locally but want to be able to deliver engaging, interactive analytic content to your stakeholders. Is there a way to do that? Lately, I’ve spent a lot of time learning and exploring tools to find the answer to that question. Luckily, the answer is yes, because there are many amazing packages that have been developed to help with this, including reactable, plotly, DT, leaflet, and crosstalk, among others. This article focuses on the latter two, and how we can build map widgets that are changeable by the user through filters, all in a self-contained HTML document.

The Use-Case: Hospital Readmissions
Constructing the Map

The Use-Case: Hospital Readmissions

A readmission occurs when a patient returns to a hospital for acute care within 30 calendar days of previously being discharged. Each year the Centers for Medicare & Medicaid Services (CMS) penalizes hospitals across the United States that have excess readmissions in a select set of high-risk patient populations. This is determined by:

Adjusting for the acuity of your hospital’s patient population by a collection of clinical risk factors
Comparing the expected readmission rate of those patients had they been treated at your hospital versus the national average hospital. The difference between those rates yields the penalty amount (with some additional steps and caveats).

For the stats gurus, it is a mixed-effects logistic regression model that uses the hospital as a random-intercept term, and clinical factors as fixed effects. You can read more about the methodology here.

The goal of this exercise is to build a map widget to identify and explore payment penalties for Wisconsin hospitals participating in the FY 2022 Hospital Readmissions Reduction Program (HRRP).

Constructing the Map

To render the map, we’ll assume you are working in an R Markdown notebook that is being rendered to HTML (as we are in this very article). The following collection of R packages will be used to do most of the heavy lifting:

Code

### Load packages
require(tidyverse)
require(leaflet)
require(crosstalk)

With those loaded, we can carry out these steps:

1. Build the dataset

We want to build a data set consisting of one row per hospital, with payment penalty information and coordinates for plotting the hospital on the map. The data needed to do this comes from a variety of sources. We can break these up into components and then put them all together into our final data set.

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:

Code

# 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.

Code

# 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)

Code

print(hrrp_results, n = 5)

# A tibble: 3,048 × 36
  `Hospital CCN` Payment adjustment f…¹ Payment Reduction Pe…² `Dual proportion`
  <chr>                           <dbl>                  <dbl>             <dbl>
1 010001                          0.990                   0.97             0.163
2 010005                          0.998                   0.15             0.141
3 010006                          0.984                   1.6              0.107
4 010007                          0.995                   0.51             0.235
5 010008                          1                       0                0.227
# ℹ 3,043 more rows
# ℹ abbreviated names: ¹`Payment adjustment factor`,
#   ²`Payment Reduction Percentage`
# ℹ 32 more variables: `Peer group assignment` <dbl>,
#   `Neutrality modifier` <dbl>, `Number of eligible discharges for AMI` <dbl>,
#   `ERR for AMI` <dbl>, `Peer group median ERR for AMI` <dbl>,
#   `Penalty indicator for AMI` <chr>, `DRG payment ratio for AMI` <dbl>, …

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:

Code

# 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
  )

Code

print(hosp_info, n = 5)

# A tibble: 5,394 × 39
  `Facility ID` `Facility Name`             Address `City/Town` State `ZIP Code`
  <chr>         <chr>                       <chr>   <chr>       <chr> <chr>     
1 010001        SOUTHEAST HEALTH MEDICAL C… 1108 R… DOTHAN      AL    36301     
2 010005        MARSHALL MEDICAL CENTERS    2505 U… BOAZ        AL    35957     
3 010006        NORTH ALABAMA MEDICAL CENT… 1701 V… FLORENCE    AL    35630     
4 010007        MIZELL MEMORIAL HOSPITAL    702 N … OPP         AL    36467     
5 010008        CRENSHAW COMMUNITY HOSPITAL 101 HO… LUVERNE     AL    36049     
# ℹ 5,389 more rows
# ℹ 33 more variables: `County/Parish` <chr>, `Telephone Number` <chr>,
#   `Hospital Type` <chr>, `Hospital Ownership` <chr>,
#   `Emergency Services` <chr>,
#   `Meets criteria for promoting interoperability of EHRs` <chr>,
#   `Meets criteria for birthing friendly designation` <chr>,
#   `Hospital overall rating` <chr>, …

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:

Code

# 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.

Code

# 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
  )

Code

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
# ℹ 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:

Code

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.

Code

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.

Code

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.

Code

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, "")
  )

Code

print(map_data, n = 5)

# A tibble: 59 × 16
  FacilityID FacilityName     Zip     lon   lat PeerGroup DualProportion Penalty
  <chr>      <chr>            <chr> <dbl> <dbl>     <dbl>          <dbl>   <dbl>
1 520002     ASPIRUS STEVENS… 54481 -89.7  44.5         4          0.262    0.12
2 520004     MAYO CLINIC HEA… 54601 -91.2  43.8         4          0.236    0.81
3 520008     WAUKESHA MEMORI… 53188 -88.2  43.0         2          0.139    0.31
4 520009     ASCENSION NE WI… 54915 -88.4  44.2         3          0.201    0.08
5 520011     MARSHFIELD MEDI… 54868 -91.7  45.5         4          0.259    0.37
# ℹ 54 more rows
# ℹ 8 more variables: AMI <chr>, COPD <chr>, HF <chr>, pneumonia <chr>,
#   CABG <chr>, `THA/TKA` <chr>, PenalizedCohortCount <dbl>,
#   PenalizedCohorts <chr>

2. Create a sharable dataset

In order to enable the map to be controlled and interacted through filters, we need to create a SharedData object using the crosstalk package. This allows the rendered HTML document to hold a connection between the data that feeds the map and the filters.

Code

# Make the shared data set
sd <- SharedData$new(data = map_data)

3. Create the base layer

The first layer of the map will consist of making a focal point on the state of Wisconsin. We will use the maps package to obtain the geometries for this region.

Code

state_outline <-
  maps::map(
    database = "state",
    regions = "wisconsin",
    fill = TRUE,
    plot = FALSE
  )

Now, we can initiate the map with the leaflet toolkit and add the first layer with a simple gray shading.

Code

my_map <- 
  leaflet() %>%
  
  # Add geographic tiles
  addTiles() %>%
  
  # Add WI state outline
  addPolygons(
    data = state_outline,
    fillColor = "gray",
    stroke = FALSE
  ) 

my_map

4. Add county lines

Next, we want to fill in the map with the Wisconsin county lines. The tigris package has a function called counties that enable us to extract these geometries.

Code

county_outlines <- 
  tigris::counties(cb = TRUE) %>%
  filter(
    STATE_NAME == "Wisconsin"
  )

We can add the county lines to the current map by supplying another layer with addPolygons. The highlightOptions argument allows us to control how a county appears when hovered over, and the label argument builds strings for tooltip display with column references to the input dataset.

Code

my_map <-
  my_map %>%
  
  # Add county outlines
  addPolygons(
    data = county_outlines,
    color = "black",
    fillColor = "#ff59f7",
    weight = 1,
    opacity = .5,
    fillOpacity = .35,
    highlightOptions = 
      highlightOptions(
        color = "black",
        weight = 3,
        bringToFront = FALSE
      ),
    label = ~NAME
  )

my_map

5. Show data points

The data points representing the hospitals can be added as yet another layer to the map. We want the color of the data points to reflect the amount of payment penalty the hospital received, so we’ll first create a color palette ranging from 0% (no penalty) to 3% (maximum penalty).

Code

# Create a color pallete
pal <- 
  colorNumeric(
    palette = "RdYlGn",
    domain = -1*seq(0,3,.1)
  )

Then we will use the shared data that was created in a prior section to add the map layer with addCircleMarkers.

Code

# Add shared data to the map
my_map %>%
  
  addCircleMarkers(
    data = sd, 
    lng = ~lon,
    lat = ~lat,
    label = ~paste0(FacilityName, " (click for info)"),
    popup = 
      ~paste0(
        "Hospital: ", FacilityName, 
        "<br>Zip Code: ", Zip,
        "<br>Peer Group: ", PeerGroup,
        "<br>Dual-Eligibility: ", round(DualProportion*100, 2), "%",
        "<br>Penalties: ", PenalizedCohorts,
        "<br>Payment Reduction: ", Penalty, "%"
      ),
    color = ~pal(-1*Penalty),
    fillOpacity = .75
  )

6. Make filters

A row of filters (and/or other objects) can be built from a call to the bscols function in the crosstalk package. This enables us to align objects in a grid across columns within rows. The first row is the widget title/header, which is just styled HTML text:

Code

# Header
bscols(
  widths = c(1, 11),
  "",
  htmltools::h3("Hospital Readmissions Reduction Program (HRRP) FY 2022 - WI Map", style = "font-style:italic;")
)

Hospital Readmissions Reduction Program (HRRP) FY 2022 - WI Map

The next row is the slider selector for filtering the penalty amount. This is built with the filter_slider function. We feed it the SharedData object previously created so that changes to the filter trigger changes to the map, which has also been supplied with the same object.

Code

# Penalty filter
bscols(
  widths = c(1, 10, 1),
  "",
  filter_slider(
    id = "penalty",
    label = "Payment Reduction",
    sharedData = sd,
    column = ~Penalty,
    step = 0.1,
    min = 0,
    max = 3,
    post = "%"
  ),
  ""
)

Payment Reduction

Similarly, the next set of rows consist of a sub-header to title the cluster of filters, and two (2) rows of dropdown selectors for each of the six (6) cohort indicators, which are analogously built with calls to filter_select:

Code

# Sub-header for cohort filters
bscols(
  widths = c(1, 11),
  "",
  "Cohort-specific penalty:"
)

Cohort-specific penalty:

Code

# Cohort filters: Row 1
bscols(
  widths = c(1, 3, 3, 3, 1),
  "",
  filter_select(
    id = "ami",
    label = "AMI",
    sharedData = sd,
    group = ~AMI
  ),
  filter_select(
    id = "cabg",
    label = "CABG",
    sharedData = sd,
    group = ~CABG
  ),
  filter_select(
    id = "copd",
    label = "COPD",
    sharedData = sd,
    group = ~COPD
  ),
  ""
)

AMI

CABG

COPD

Code

# Cohort filters: Row 2
bscols(
  widths = c(1, 3, 3, 3, 1),
  "",
  filter_select(
    id = "hf",
    label = "HF",
    sharedData = sd,
    group = ~HF
  ),
  filter_select(
    id = "pneu",
    label = "PNEU",
    sharedData = sd,
    group = ~pneumonia
  ),
  filter_select(
    id = "tja",
    label = "THA/TKA",
    sharedData = sd,
    group = ~`THA/TKA`
  ),
  ""
)

PNEU

THA/TKA

7. Format order of appearance

Notice that when you adjust the filters in the last section, they change the map from step 5. This is because they are all referencing the same SharedData object. So, our final step is that we simply need to rearrange the execution order of the code so that the objects appear in the relative positions we would like. The code chunk creating the following output is a copy-paste of what we executed in steps 2-6, just with some rearranging. We will set echo = FALSE in the code chunk options so that only the output renders together. Note that we created a different SharedData object here so that the filters in the prior steps don’t impact this output.

Hospital Readmissions Reduction Program (HRRP) FY 2022 - WI Map

Payment Reduction

Cohort-specific penalty:

AMI

CABG

COPD

PNEU

THA/TKA

--- title: "Building a filterable map with {leaflet} and {crosstalk} for hospital readmissions" description: "For runtime-free analytical content delivery" author: "Alex Zajichek" date: "6/18/2022" image: "feature.png" categories: - Healthcare - Maps - Javascript format: html: code-fold: true code-tools: true --- Want to make this map widget? Read On! `r emo::ji("backhand index pointing down")``r emo::ji("backhand index pointing down")``r emo::ji("backhand index pointing down")` ```{r, echo = FALSE, message = FALSE, warning = FALSE, error = FALSE, comment = FALSE, results = 'hide'} ### Load packages require(tidyverse) require(leaflet) require(crosstalk) ### 1. Import CMS data ## i. FY22 final rule (indicates penalties for each hospital) # 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) # 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) ## ii. Hospital Information # Set path to file (from CMS) hosp_file <- "/Users/alexzajichek/Documents/GitHub/Hospital_General_Information.csv" # Import the file hosp_info <- read_csv( file = hosp_file ) ### 2. Obtain zip code centroids # Get set of zip codes for Wisconsin zips_wi <- tigris::zctas(year = 2010, state = "WI") # 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 ) ### 3. Build an analysis dataset 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` ) %>% # 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 ) ) %>% # 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" ) %>% # 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, "") ) ### 4. Map data components # State outline state_outline <- maps::map( database = "state", regions = "wisconsin", fill = TRUE, plot = FALSE ) # County outlines county_outlines <- tigris::counties(cb = TRUE) %>% filter( STATE_NAME == "Wisconsin" ) ### 4. Build map ## i. State outline my_map <- leaflet() %>% # Add geographic tiles addTiles() %>% # Add WI state outline addPolygons( data = state_outline, fillColor = "gray", stroke = FALSE ) ## ii. County outline my_map <- my_map %>% # Add county outlines addPolygons( data = county_outlines, color = "black", fillColor = "#ff59f7", weight = 1, opacity = .5, fillOpacity = .35, highlightOptions = highlightOptions( color = "black", weight = 3, bringToFront = FALSE ), label = ~NAME ) ``` ```{r, echo = FALSE, message = FALSE, warning = FALSE, error = FALSE, comment = FALSE} ## iii. Add hospitals # Make the shared data set sd <- SharedData$new(data = map_data) sd2 <- SharedData$new(data = map_data) # Make filters bscols( widths = c(1, 11), "", htmltools::h3("Hospital Readmissions Reduction Program (HRRP) FY 2022 - WI Map", style = "font-style:italic;") ) # Row 1: bscols( widths = c(1, 10, 1), "", filter_slider( id = "penalty2", label = "Payment Reduction", sharedData = sd2, column = ~Penalty, step = 0.1, min = 0, max = 3, post = "%" ), "" ) # Row 2-3: bscols( widths = c(1, 11), "", "Cohort-specific penalty:" ) bscols( widths = c(1, 3, 3, 3, 1), "", filter_select( id = "ami2", label = "AMI", sharedData = sd2, group = ~AMI ), filter_select( id = "cabg2", label = "CABG", sharedData = sd2, group = ~CABG ), filter_select( id = "copd2", label = "COPD", sharedData = sd2, group = ~COPD ), "" ) bscols( widths = c(1, 3, 3, 3, 1), "", filter_select( id = "hf2", label = "HF", sharedData = sd2, group = ~HF ), filter_select( id = "pneu2", label = "PNEU", sharedData = sd2, group = ~pneumonia ), filter_select( id = "tja2", label = "THA/TKA", sharedData = sd2, group = ~`THA/TKA` ), "" ) # Create a color pallete pal <- colorNumeric( palette = "RdYlGn", domain = -1*seq(0,3,.1) ) # Add shared data to the map my_map %>% addCircleMarkers( data = sd2, # A dataset containing coordinates for each zip code (or anything) lng = ~lon, # The coordinates of that entitity lat = ~lat, label = ~paste0(FacilityName, " (click for info)"), popup = ~paste0( "Hospital: ", FacilityName, " Zip Code: ", Zip, " Peer Group: ", PeerGroup, " Dual-Eligibility: ", round(DualProportion*100, 2), "%", " Penalties: ", PenalizedCohorts, " Payment Reduction: ", Penalty, "%" ), color = ~pal(-1*Penalty), fillOpacity = .75 ) ``` So you run `R` locally but want to be able to deliver engaging, interactive analytic content to your stakeholders. Is there a way to do that? Lately, I've spent a lot of time learning and exploring tools to find the answer to that question. Luckily, the answer is _yes_, because there are many amazing packages that have been developed to help with this, including [`reactable`](https://glin.github.io/reactable/), [`plotly`](https://plotly.com/r/), [`DT`](https://rstudio.github.io/DT/), [`leaflet`](https://rstudio.github.io/leaflet/), and [`crosstalk`](https://rstudio.github.io/crosstalk/), among others. This article focuses on the latter two, and how we can build map widgets that are changeable by the user through filters, all in a self-contained HTML document. # Table of Contents * [The Use-Case: Hospital Readmissions](#usecase) * [Constructing the Map](#mapoverview) 1. [Build the dataset](#usecasedataset) 2. [Create a shareable dataset](#shareddata) 3. [Create the base layer](#baselayer) 4. [Add county lines](#countylines) 5. [Show data points](#datapoints) 6. [Make filters](#filters) 7. [Format order of appearance](#formatorder) # The Use-Case: Hospital Readmissions {#usecase} A _readmission_ occurs when a patient returns to a hospital for acute care within 30 calendar days of previously being discharged. Each year the [_Centers for Medicare & Medicaid Services (CMS)_](https://www.cms.gov/) [penalizes hospitals](https://www.cms.gov/Medicare/Medicare-Fee-for-Service-Payment/AcuteInpatientPPS/Readmissions-Reduction-Program) across the United States that have _excess_ readmissions in a select set of high-risk patient populations. This is determined by: 1. Adjusting for the acuity of your hospital's patient population by a collection of clinical risk factors 2. Comparing the expected readmission rate of _those_ patients had they been treated at _your_ hospital versus the _national average_ hospital. The difference between those rates yields the penalty amount (with some additional steps and caveats). For the stats gurus, it is a mixed-effects logistic regression model that uses the hospital as a random-intercept term, and clinical factors as fixed effects. You can read more about the methodology [here](https://qualitynet.cms.gov/inpatient/measures/readmission/methodology). The goal of this exercise is to build a map widget to identify and explore payment penalties for Wisconsin hospitals participating in the _FY 2022 Hospital Readmissions Reduction Program (HRRP)_. # Constructing the Map {#mapoverview} To render the map, we'll assume you are working in an [R Markdown](https://rmarkdown.rstudio.com/) notebook that is being rendered to HTML (as we are in this very article). The following collection of `R` packages will be used to do most of the heavy lifting: ```{r, eval = FALSE} ### Load packages require(tidyverse) require(leaflet) require(crosstalk) ``` With those loaded, we can carry out these steps: ## 1. Build the dataset {#usecasedataset} We want to build a data set consisting of one row per hospital, with payment penalty information and coordinates for plotting the hospital on the map. The data needed to do this comes from a variety of sources. We can break these up into components and then put them all together into our final data set. ### i. Payment penalties The penalty amounts for each hospital are documented in the [CMS FY2022 Final Rule](https://www.cms.gov/medicare/acute-inpatient-pps/fy-2022-ipps-final-rule-home-page). The file we need is contained within a `zip` file; we can use [this](rpubs.com/otienodominic/398952) solution for motivation of our approach. First, we can download the file: ```{r, eval = FALSE} # 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. ```{r, eval = FALSE} # 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) ``` ```{r} print(hrrp_results, n = 5) ``` ### 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](https://data.cms.gov/provider-data/dataset/xubh-q36u), which we can import from directly: ```{r, eval = FALSE} # 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 ) ``` ```{r} print(hosp_info, n = 5) ``` ### 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: ```{r, eval = FALSE} # 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. ```{r, eval = FALSE} # 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 ) ``` ```{r} print(zips_wi_centroids, n = 5) ``` ### 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: ```{r, eval = FALSE} 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. ```{r, eval = FALSE} 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. ```{r, eval = FALSE} 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. ```{r, eval = FALSE} 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, "") ) ``` ```{r} print(map_data, n = 5) ``` ## 2. Create a sharable dataset {#shareddata} In order to enable the map to be controlled and interacted through filters, we need to create a `SharedData` object using the `crosstalk` package. This allows the rendered HTML document to hold a connection between the data that feeds the map and the filters. ```{r, eval = FALSE} # Make the shared data set sd <- SharedData$new(data = map_data) ``` ## 3. Create the base layer {#baselayer} The first layer of the map will consist of making a focal point on the state of Wisconsin. We will use the `maps` package to obtain the geometries for this region. ```{r, eval = FALSE} state_outline <- maps::map( database = "state", regions = "wisconsin", fill = TRUE, plot = FALSE ) ``` Now, we can initiate the map with the `leaflet` toolkit and add the first layer with a simple gray shading. ```{r} my_map <- leaflet() %>% # Add geographic tiles addTiles() %>% # Add WI state outline addPolygons( data = state_outline, fillColor = "gray", stroke = FALSE ) my_map ``` ## 4. Add county lines {#countylines} Next, we want to fill in the map with the Wisconsin county lines. The `tigris` package has a function called `counties` that enable us to extract these geometries. ```{r, eval = FALSE} county_outlines <- tigris::counties(cb = TRUE) %>% filter( STATE_NAME == "Wisconsin" ) ``` We can add the county lines to the current map by supplying another layer with `addPolygons`. The `highlightOptions` argument allows us to control how a county appears when hovered over, and the `label` argument builds strings for tooltip display with column references to the input dataset. ```{r, warning = FALSE} my_map <- my_map %>% # Add county outlines addPolygons( data = county_outlines, color = "black", fillColor = "#ff59f7", weight = 1, opacity = .5, fillOpacity = .35, highlightOptions = highlightOptions( color = "black", weight = 3, bringToFront = FALSE ), label = ~NAME ) my_map ``` ## 5. Show data points {#datapoints} The data points representing the hospitals can be added as yet another layer to the map. We want the color of the data points to reflect the amount of payment penalty the hospital received, so we'll first create a color palette ranging from 0% (no penalty) to 3% (maximum penalty). ```{r, eval = FALSE} # Create a color pallete pal <- colorNumeric( palette = "RdYlGn", domain = -1*seq(0,3,.1) ) ``` Then we will use the shared data that was created in a [prior section](#shareddata) to add the map layer with `addCircleMarkers`. ```{r} # Add shared data to the map my_map %>% addCircleMarkers( data = sd, lng = ~lon, lat = ~lat, label = ~paste0(FacilityName, " (click for info)"), popup = ~paste0( "Hospital: ", FacilityName, " Zip Code: ", Zip, " Peer Group: ", PeerGroup, " Dual-Eligibility: ", round(DualProportion*100, 2), "%", " Penalties: ", PenalizedCohorts, " Payment Reduction: ", Penalty, "%" ), color = ~pal(-1*Penalty), fillOpacity = .75 ) ``` ## 6. Make filters {#filters} A row of filters (and/or other objects) can be built from a call to the `bscols` function in the `crosstalk` package. This enables us to align objects in a grid across columns within rows. The first row is the widget title/header, which is just styled HTML text: ```{r} # Header bscols( widths = c(1, 11), "", htmltools::h3("Hospital Readmissions Reduction Program (HRRP) FY 2022 - WI Map", style = "font-style:italic;") ) ``` The next row is the slider selector for filtering the penalty amount. This is built with the `filter_slider` function. We feed it the `SharedData` object [previously created](#shareddata) so that changes to the filter trigger changes to the map, which has also been supplied with the same object. ```{r} # Penalty filter bscols( widths = c(1, 10, 1), "", filter_slider( id = "penalty", label = "Payment Reduction", sharedData = sd, column = ~Penalty, step = 0.1, min = 0, max = 3, post = "%" ), "" ) ``` Similarly, the next set of rows consist of a sub-header to title the cluster of filters, and two (2) rows of dropdown selectors for each of the six (6) cohort indicators, which are analogously built with calls to `filter_select`: ```{r} # Sub-header for cohort filters bscols( widths = c(1, 11), "", "Cohort-specific penalty:" ) # Cohort filters: Row 1 bscols( widths = c(1, 3, 3, 3, 1), "", filter_select( id = "ami", label = "AMI", sharedData = sd, group = ~AMI ), filter_select( id = "cabg", label = "CABG", sharedData = sd, group = ~CABG ), filter_select( id = "copd", label = "COPD", sharedData = sd, group = ~COPD ), "" ) # Cohort filters: Row 2 bscols( widths = c(1, 3, 3, 3, 1), "", filter_select( id = "hf", label = "HF", sharedData = sd, group = ~HF ), filter_select( id = "pneu", label = "PNEU", sharedData = sd, group = ~pneumonia ), filter_select( id = "tja", label = "THA/TKA", sharedData = sd, group = ~`THA/TKA` ), "" ) ``` ## 7. Format order of appearance {#formatorder} Notice that when you adjust the filters in the [last section](#filters), they change the map from [step 5](#datapoints). This is because they are all referencing the same `SharedData` object. So, our final step is that we simply need to rearrange the execution order of the code so that the objects appear in the relative positions we would like. The code chunk creating the following output is a copy-paste of what we executed in steps 2-6, just with some rearranging. We will set `echo = FALSE` in the code chunk options so that only the output renders together. _Note that we created a different `SharedData` object here so that the filters in the prior steps don't impact this output._ ```{r, echo = FALSE} # Make the shared data set sd3 <- SharedData$new(data = map_data) # Header bscols( widths = c(1, 11), "", htmltools::h3("Hospital Readmissions Reduction Program (HRRP) FY 2022 - WI Map", style = "font-style:italic;") ) # Penalty filter bscols( widths = c(1, 10, 1), "", filter_slider( id = "penalty", label = "Payment Reduction", sharedData = sd3, column = ~Penalty, step = 0.1, min = 0, max = 3, post = "%" ), "" ) # Sub-header for cohort filters bscols( widths = c(1, 11), "", "Cohort-specific penalty:" ) # Cohort filters: Row 1 bscols( widths = c(1, 3, 3, 3, 1), "", filter_select( id = "ami", label = "AMI", sharedData = sd3, group = ~AMI ), filter_select( id = "cabg", label = "CABG", sharedData = sd3, group = ~CABG ), filter_select( id = "copd", label = "COPD", sharedData = sd3, group = ~COPD ), "" ) # Cohort filters: Row 2 bscols( widths = c(1, 3, 3, 3, 1), "", filter_select( id = "hf", label = "HF", sharedData = sd3, group = ~HF ), filter_select( id = "pneu", label = "PNEU", sharedData = sd3, group = ~pneumonia ), filter_select( id = "tja", label = "THA/TKA", sharedData = sd3, group = ~`THA/TKA` ), "" ) # Add shared data to the map my_map %>% addCircleMarkers( data = sd3, lng = ~lon, lat = ~lat, label = ~paste0(FacilityName, " (click for info)"), popup = ~paste0( "Hospital: ", FacilityName, " Zip Code: ", Zip, " Peer Group: ", PeerGroup, " Dual-Eligibility: ", round(DualProportion*100, 2), "%", " Penalties: ", PenalizedCohorts, " Payment Reduction: ", Penalty, "%" ), color = ~pal(-1*Penalty), fillOpacity = .75 ) ```