## data
ra_ft = structure(list(jurisdiction_id = c(258L, 258L, 258L, 258L, 258L, 258L, 258L, 258L, 258L, 258L, 258L, 258L, 258L, 258L, 258L), hazard_name = c("Pandemic Influenza", "Hurricane/Tropical Storm", "Extreme Heat", "Nuclear Attack", "Drought ", "Flood", "Windstorm ", "Biological Disease Outbreak", "Power Failure ", "Biological Terrorism - Communicable (including A - B - C agents)", "Mass Casualty Incidents", "Biological Terrorism - Non-Communicable (including A - B - C agents)", "Water Supply Contamination - environmental", "Food Borne Disease", "Extreme Cold"), hazard_risk_index = c(152.17, 139.72, 82.09, 55.84, 56.81, 61.47, 58.19, 60.4, 51.72, 49.37, 53.08, 45.57, 47.89, 50.29, 47.83), residual_risk_index = c(11.22, 9.62, 6.08, 5.1, 4.76, 4.53, 4.52, 4.41, 4.41, 4.08, 3.93, 3.82, 3.82, 3.59, 3.49), probability_score = c(3, 3.22, 2.65, 1, 2.5, 1.87, 2.62, 1.33, 2.29, 1, 1.72, 1, 1.18, 1.3, 2.53)), row.names = c(NA, -15L), class = c("tbl_df", "tbl", "data.frame"))

library(tidyverse)
prob_factor = with(ra_ft, max(hazard_risk_index) / max(probability_score))

long_dat = ra_ft |>
  mutate(hazard_name = reorder(hazard_name, -residual_risk_index)) |>
  pivot_longer(cols = c("hazard_risk_index", "residual_risk_index", "probability_score"))

ggplot(
    data = long_dat |> filter(name != "probability_score"),
    aes(x = hazard_name, y = value, fill = name)
  ) +
  geom_col(
    #aes(fill = name),
    position = position_dodge()
  ) +
  geom_point(
    data = long_dat |> filter(name == "probability_score"),
    aes(y = value * prob_factor),
    #color = "red4"
    shape = 21
  ) +
  scale_fill_manual(
    values = c(
      hazard_risk_index = "cadetblue2",
      residual_risk_index = "yellow3",
      probability_score = "red4"
    ),
    limits = c(
      "hazard_risk_index",
      "residual_risk_index",
      "probability_score"
    ),
    labels = \(x) {x |>  str_replace_all("_", " ") |> str_to_title()}
  ) +
  scale_y_continuous(
    name = "Risk Indices",
    sec.axis = sec_axis(~ . / prob_factor, name = "Probability Score"),
    expand = expansion(mult = c(0, 0.1))
  ) +
  scale_x_discrete(
    labels = \(x) str_replace(x, pattern = fixed("("), replacement = "\n(")
  ) +
  labs(
    x = NULL,
    fill = NULL,
    title = "Hazard Probability, Residual Risk, and Probability Score"
  ) +
  ## putting the legend at the bottom this way requires
  ## ggplot2 version 3.5.0 or higher
  guides(fill = guide_legend(
    position = "bottom",
    label.position = "bottom"
  )) +
  theme_bw() +
  theme(
    panel.grid.major.x = element_blank(),
    panel.border = element_blank(),
    axis.text.x = element_text(angle = 45, hjust = 1),
    axis.ticks = element_line(color = "gray85"),
    plot.margin = margin(t = 5.5, r = 5.5, b = 10, l = 20, unit = "pt")
  )