Functions

The code in this script performs a number of useful tasks required for the analysis.

Source code

The full R script is available at:

• R/functions.R

This R script file is itself reliant on the following helper files:

• R/setup.R

Drawing 60 minute isochrones around each transplant center

The following function uses the mapboxapi package to calculate 60-minute travel isochrones around each transplant center.

The function assumes the that you have an API token configured in R. For details of how to do this, see the package documentation.

Its inputs are as follows:

• SFObject: This is an object of class sf that contains the location of each transplant center
• year: year as a character. Must be a year in yearlist defined in R/setup.R.

#Define a function that takes an SF data object and returns 60 minute isochrones
set_60_min_isochrone_at_7AM<-function(SFObject, year) {
  
  if (nrow(SFObject) == 0) {
    return(
      SFObject %>%
        st_transform(5070)
    )
  }else{
    tempDF<-SFObject%>%
      mutate(TimeZone=tz_lookup(., method="accurate"))%>%
      st_transform(4326)
    
    FinalDF <- vector("list", length(timezones))
    
    
    for(i in 1:length(timezones))
    {
      TZ<-converted_times[[year]][[i,1]]
      
      UTCTime<-format(converted_times[[year]][[i,3]], "%Y-%m-%dT%H:%M:%SZ")
      
      if (nrow(filter(tempDF,TimeZone==TZ))>0)
      {
        temp<-tempDF%>%
          filter(TimeZone==TZ)%>%
          mb_isochrone(time = 60, 
                       profile = "driving-traffic",
                       depart_at = UTCTime,
                       id_column = "OTCCode")%>%
          rename(OTCCode=id)
        FinalDF[[i]]<-temp
      }
    }
    
    FinalDF%>%
      bind_rows()%>%
      st_transform(5070)%>%
      return()
  }
}

Converting data frame of summary results into a gt format table

This function takes a summary results data frame and formats it into a publication-ready table using the gt package. The input data are assumed to contain one row per characteristic and year, with separate columns for the numerator, denominator, and percentage of the population meeting the specified access criterion.

The function first reshapes the data from long to wide format so that each study year appears as a grouped set of columns. It then converts the Characteristic column into a factor with a user-specified ordering and labeling, ensuring that rows appear in a meaningful and reproducible sequence rather than default alphabetical order.

After restructuring the data, the function applies a consistent gt table layout. Columns are grouped under year-specific spanner headers, counts are displayed with thousands separators, and percentages are rounded to one decimal place. The table title is supplied as an argument, allowing the same function to be reused across multiple result tables.

To improve readability, the function also applies subtle background shading to distinguish columns belonging to different years and adds vertical borders between year groupings. The result is a standardized comparative table that can be reused throughout the project for presenting access estimates across organs, center classifications, and distance definitions.

This function also makes repeated use of purrr::reduce(), which is a functional programming tool for applying the same transformation cumulatively across a sequence of values. Here, the sequence is the set of study years. Each iteration takes the current gt table and adds another layer of formatting, such as a year-specific spanner, fill color, or border. This avoids repetitive code and makes the table-building process both more compact and more generalizable to analyses involving different numbers of years.

convert_resultsdf_to_table<-function(resultsdf, 
                                     year_list_fn=year_list,
                                     row_values,
                                     row_labels,
                                     table_title){
  
  
  resultsdf%>%
    
    #Reshape data set
    pivot_wider(
      names_from = Year,
      values_from = c(Numerator, Denominator, Percentage),
      names_vary = "slowest"
    )%>%
    
    #Use the table defined above to turn Characteristic into a factor variable
    mutate(Characteristic=factor(Characteristic,
                                 levels = row_values,
                                 labels = row_labels))%>%
    arrange(Characteristic)%>%
    
    #Convert to a gt format table
    gt()%>%
    
    #Groups columns by year with a spanner label
    {
      purrr::reduce(
        year_list_fn,
        ~ .x %>%
          tab_spanner(
            label = .y,
            columns = matches(paste0("_", .y, "$"))
          ),
        .init = .
      )
    } %>%
    
    #Sets coluns to round to 0 decimal places for numbers and 1 decimal place for percentages
    fmt_number(
      columns = matches(paste0("Numerator_|Denominator_")),
      decimals = 0,
      use_seps = TRUE
    ) %>%
    fmt_number(
      columns = matches("Percentage_"),
      decimals = 1
    )%>%tab_header(
      title = table_title
    )%>%
    
    #Removes "_year" from columns
    cols_label_with(~ "Numerator",   columns = matches("^Numerator_")) %>%
    cols_label_with(~ "Denominator", columns = matches("^Denominator_")) %>%
    cols_label_with(~ "Percentage",  columns = matches("^Percentage_"))%>%
    
    #Shades columns by year for easier reading
    
    {
      purrr::reduce(
        seq_along(year_list_fn),
        ~ .x %>%
          tab_style(
            style = cell_fill(
              color = colorRampPalette(c("lightgray","lightblue"))(
                length(year_list_fn)
              )[.y]            ),
            locations = cells_body(
              columns = matches(paste0("_", year_list_fn[.y], "$"))
            )
          ),
        .init = .
      )
    }%>%
    
    #Add separators between years
    {
      purrr::reduce(
        seq_len(length(year_list_fn) - 1),
        ~ .x %>%
          tab_style(
            style = cell_borders(
              sides = "right",
              weight = px(2),
              color = "black"
            ),
            locations = cells_body(
              columns = matches(paste0("Percentage_", year_list_fn[.y], "$"))
            )
          ) %>%
          tab_style(
            style = cell_borders(
              sides = "left",
              weight = px(2),
              color = "black"
            ),
            locations = cells_body(
              columns = matches(paste0("Numerator_", year_list_fn[.y + 1], "$"))
            )
          ),
        .init = .
      )
    }
  
  
}


}

Creates paired plots of maps with percentages in the plot

This function creates side-by-side map panels comparing geographic access in two different years, while also embedding the corresponding summary percentages directly into the plot titles. For a given organ, distance definition, and center category, the function pulls the relevant summary statistics from the precomputed results tables and combines them with the geographic buffer objects used to draw the maps. This allows each figure to show both the spatial distribution of access and the headline percentage values in a single visual.

The function is designed to compare two years, by default 2017 and 2022. For each year, it overlays transplant center catchment areas in red on a national base map, while blue dots represent the estimated distribution of people living with HIV. The title above each panel reports the percentage, numerator, and denominator for the selected outcome, allowing the reader to interpret the geographic pattern alongside the underlying quantitative summary.

An optional argument, include_nonHIV, allows the figure title to also display the corresponding results for people without HIV. This makes it possible to present HIV and non-HIV comparisons within the same figure when desired, while keeping the map itself visually simple. The function also performs basic checks to ensure that the requested years and outcome labels are valid and that the buffer object supplied is likely to match the intended analysis.

Finally, the two maps are combined into a single annotated figure using cowplot, with a shared title and a brief explanatory note. The result is a presentation-style comparative graphic that is useful for showing changes over time in transplant center catchment areas and in the proportion of people with HIV who live within those areas.

make_paired_plot<-function(organ,
                           distance,
                           outcome,
                           buffer_list,
                           year1="2017",
                           year2="2022",
                           plottitle,
                           include_nonHIV=FALSE)
{
  
  # ----Error checking----
  
  # Validate include_nonHIV
  if (!is.logical(include_nonHIV) || length(include_nonHIV) != 1 || is.na(include_nonHIV)) {
    stop("include_nonHIV must be a single TRUE or FALSE value.")
  }
  
  
  #Access name of object passed to `buffer_list`
  buffer_name <- deparse(substitute(buffer_list))
  
  #Make sure name of outcome matches buffer_list
  if (!stringr::str_detect(
    stringr::str_to_lower(buffer_name),
    stringr::str_to_lower(outcome)
  )) {
    warning(glue::glue(
      "Buffer object ('{buffer_name}') may not match outcome '{outcome}'."
    ))
  }
  
  if (!year1 %in% year_list) {
    stop(glue::glue("year1 ({year1}) is not in year_list."))
  }
  
  if (!year2 %in% year_list) {
    stop(glue::glue("year2 ({year2}) is not in year_list."))
  }
  
  if (!outcome %in% c("active", "HIV", "HOPE")) {
    stop(glue::glue("outcome ({outcome}) is not in c('active', 'HIV', 'HOPE')."))
  }
  
  # ----Main function----
  
  #Keep only elements of Results_HIV_df necessary for plot
  summary_df_HIV<-Results_HIV_df[[organ]]%>%
    filter(str_detect(Characteristic, outcome))%>%
    filter(str_detect(Characteristic, distance))
  
  #parameters of plot assigned
  parameter_plot<-list()
  parameter_plot$y_title<-0.9
  parameter_plot$y_footnote<-0.1
  
  if (include_nonHIV)
  {
    parameter_plot$y_title<-0.97
    parameter_plot$y_footnote<-0.01
  }
  
  
  #Keep only elements of Results_nonHIV_df necessary for plot, if even needed in the first place
  if (include_nonHIV)
  {
    summary_df_nonHIV<-Results_nonHIV_df[[organ]]%>%
      filter(str_detect(Characteristic, outcome))%>%
      filter(str_detect(Characteristic, distance))
  }
  
  #Creates text for nonHIV (is NULL if )
  title_text_nonHIV<-list()
  
  if (include_nonHIV==TRUE) {
    title_text_nonHIV[[year1]]<-
      glue("<br>",
           "<span style='font-size:30pt; font-weight:normal;'>HIV negative:</span> ",
           "<span style='font-size:30pt; font-weight:bold;'>",
           "{summary_df_nonHIV$Percentage[summary_df_nonHIV$Year == as.numeric(year1)]}%",
           "</span>",
           "<br>",
           "{comma(summary_df_nonHIV$Numerator[summary_df_nonHIV$Year == as.numeric(year1)])}",
           " / ",
           "{comma(summary_df_nonHIV$Denominator[summary_df_nonHIV$Year == as.numeric(year1)])}")
    
    title_text_nonHIV[[year2]]<-
      glue("<br>",
           "<span style='font-size:30pt; font-weight:normal;'>HIV negative:</span> ",
           "<span style='font-size:30pt; font-weight:bold;'>",
           "{summary_df_nonHIV$Percentage[summary_df_nonHIV$Year == as.numeric(year2)]}%",
           "</span>",
           "<br>",
           "{comma(summary_df_nonHIV$Numerator[summary_df_nonHIV$Year == as.numeric(year2)])}",
           " / ",
           "{comma(summary_df_nonHIV$Denominator[summary_df_nonHIV$Year == as.numeric(year2)])}")
  }
  else {
    title_text_nonHIV[[year1]]<-""
    title_text_nonHIV[[year2]]<-""
  }
  
  
  
  #Left graph
  Plot1<-ggplot() +
    geom_sf(data = county_dots_transformed[[year1]], size = 0.3, color = "navy", alpha = 0.5) +
    theme_void() +
    theme(
      axis.title.y = element_text(size = 36, angle = 90, vjust = 0.5),
      plot.title = element_textbox(size = 16, hjust = 0.5,
                                   halign = 0.5)
    )+
    labs(y=year1,
         title = glue("<span style='font-size:30pt; font-weight:normal;'>PLWH:</span> ",
                      "<span style='font-size:30pt; font-weight:bold;'>",
                      "{summary_df_HIV$Percentage[summary_df_HIV$Year == as.numeric(year1)]}%",
                      "</span>",
                      "<br>",
                      "{comma(summary_df_HIV$Numerator[summary_df_HIV$Year == as.numeric(year1)])}",
                      " / ",
                      "{comma(summary_df_HIV$Denominator[summary_df_HIV$Year == as.numeric(year1)])}",
                      title_text_nonHIV[[year1]]
         ))+
    geom_sf(data=buffer_list[[organ]][[distance]][[year1]], color="red", fill=NA)+
    geom_sf(data=states_sf_transformed, fill=NA)
  
  
  #Right plot
  Plot2<-ggplot() +
    geom_sf(data = county_dots_transformed[[year2]], size = 0.3, color = "navy", alpha = 0.5) +
    theme_void() +
    theme(
      axis.title.y = element_text(size = 36, angle = 90, vjust = 0.5),
      plot.title = element_textbox(size = 16, hjust = 0.5,
                                   halign = 0.5)
    )+
    labs(y=year2,
         title = glue("<span style='font-size:30pt; font-weight:normal;'>PLWH:</span> ",
                      "<span style='font-size:30pt; font-weight:bold;'>",
                      "{summary_df_HIV$Percentage[summary_df_HIV$Year == as.numeric(year2)]}%",
                      "</span>",
                      "<br>",
                      "{comma(summary_df_HIV$Numerator[summary_df_HIV$Year == as.numeric(year2)])}",
                      " / ",
                      "{comma(summary_df_HIV$Denominator[summary_df_HIV$Year == as.numeric(year2)])}",
                      title_text_nonHIV[[year2]]
         ))+
    geom_sf(data=buffer_list[[organ]][[distance]][[year2]], color="red", fill=NA)+
    geom_sf(data=states_sf_transformed, fill=NA)
  
  #Combine plots
  combined<-cowplot::plot_grid(Plot1, Plot2,
                               ncol = 2, align = "hv")
  
  final_plot <- ggdraw() +
    draw_plot(combined) +
    draw_label(
      plottitle,
      x = 0.5, y = parameter_plot$y_title,   # top center
      hjust = 0.5, vjust = 1,   # anchor to the top edge
      fontface = 'bold',
      size = 25
    )+
    # 
    draw_label(
      "Each blue dot represents 100 people with HIV\nRed outlines represent the catchment area for each transplant center",
      x = 0.02, y = parameter_plot$y_footnote,          # left margin position
      hjust = 0, vjust = 0,        # align text to the left edge
      size = 14,
      fontface = "italic"
    )
  
  return(final_plot)
  
}

Plot the number of people with HIV in the catchment area of each transplant center

This function uses the ggbeeswarm package to make GraphPad Prism-type plots comparing the number of people living with HIV in the catchment area of each transplant center in the years specified.

plot_center_HIV_catchment<-function(organ,
                                    distance,
                                    buffer_list,
                                    year1="2017",
                                    year2="2022"){
  
  
  
  # ----Error checking----
  
  #Access name of object passed to `buffer_list`
  buffer_name <- deparse(substitute(buffer_list))
  
  if (!year1 %in% year_list) {
    stop(glue::glue("year1 ({year1}) is not in year_list."))
  }
  
  if (!year2 %in% year_list) {
    stop(glue::glue("year2 ({year2}) is not in year_list."))
  }
  
  
  
  temp_DF<-list()
  
  temp_DF[[year1]]<-buffer_list%>%
    pluck(year1)%>%
    st_drop_geometry()%>%
    group_by(OTCCode)%>%
    summarize(total_cases=round(sum(tract_cases), 0))%>%
    mutate(year=as.numeric(year1))%>%
    mutate(Center_classification=case_when(
      OTCCode %in% HOPE_center_volumes[[organ]][[year1]]$REC_CTR_CD~"HOPE center",
      OTCCode %in% HIV__center_volumes[[organ]][[year1]]$REC_CTR_CD~"≥1 HIV R+ transplants",
      OTCCode %in% Active_tx_centers[[organ]][[year1]]~"No HIV R+ transplants"
    ))%>%
    filter(!is.na(OTCCode))
  
  
  temp_DF[[year2]]<-buffer_list%>%
    pluck(year2)%>%
    st_drop_geometry()%>%
    group_by(OTCCode)%>%
    summarize(total_cases=round(sum(tract_cases), 0))%>%
    mutate(year=as.numeric(year2))%>%
    mutate(Center_classification=case_when(
      OTCCode %in% HOPE_center_volumes[[organ]][[year2]]$REC_CTR_CD~"HOPE center",
      OTCCode %in% HIV__center_volumes[[organ]][[year2]]$REC_CTR_CD~"≥1 HIV R+ transplants",
      OTCCode %in% Active_tx_centers[[organ]][[year2]]~"No HIV R+ transplants"
    ))%>%
    filter(!is.na(OTCCode))
  
  combined_DF<-bind_rows(temp_DF)%>%
    mutate(Center_classification=factor(Center_classification,
                                        levels=c("HOPE center","≥1 HIV R+ transplants", "No HIV R+ transplants"),
                                        labels=c("HOPE center","≥1 HIV R+ transplants", "No HIV R+ transplants")
    ))%>%
    mutate(Center_classification_color=case_when(
      OTCCode %in% HOPE_center_volumes[[organ]][[year1]]$REC_CTR_CD~"HOPE center",
      OTCCode %in% HIV__center_volumes[[organ]][[year1]]$REC_CTR_CD~"≥1 HIV R+ transplants",
      OTCCode %in% Active_tx_centers[[organ]][[year1]]~"No HIV R+ transplants",
      TRUE~glue("No transplants in {year1}")
    ))
  
  combined_DF %>% ggplot()+
    geom_violin(mapping = aes(x=Center_classification, y=total_cases), color = NA, alpha=0.5,fill = "grey", 
                draw_quantiles = 1)+
    stat_summary(mapping = aes(x=Center_classification, y=total_cases),fun = median, 
                 fun.min = median, fun.max = median, 
                 geom = "crossbar", width = 0.9, fatten=1.5)+
    geom_quasirandom(mapping = aes(x=Center_classification, y=total_cases, color=Center_classification_color), 
                     cex=2, alpha=0.7)+
    stat_compare_means(method = "kruskal.test",
                       mapping = aes(x = Center_classification, y = total_cases)) +
    scale_y_continuous(labels = scales::comma)+
    guides(color = guide_legend(title.position = "top", nrow = 1))+
    theme_classic()+
    theme(legend.position = "bottom",
          legend.title = element_text(hjust = 0.5),  # centers title
          legend.box = "horizontal")+
    theme(axis.text.x = element_text(
      angle = 45,
      hjust = 1,
      vjust = 1))+
    facet_grid(.~year)+
    labs(x = "Center HIV R+ volume", y="PLWH in catchment area", color=glue("{year1} status"))+
    scale_color_discrete(breaks = c("HOPE center",
                                    "≥1 HIV R+ transplants", 
                                    "No HIV R+ transplants",
                                    glue("No transplants in {year1}"))
    )
  
  
  
}