# pertussis_analysis.R
#
# Reference code for automating an Excel-based analysis of an MDSS export of a
# pertussis query result. Written for MI-Support's 24 July 2026 technical
# workshop in Lansing, MI
#
# Authors: Jack Jacobs (jckjcbs@umich.edu)
#          MI-Support (micom-misupport@umich.edu)


#### User inputs ####

# Only change the values in this section

# Location of the MDSS export file
mdss_path <- "data/FAKE_pertussis_mdss.csv"

# Location of the Ingham County annual population by sex file
population_path <- "data/population.csv"

# Location where you want to export the cases per month by sex table
month_sex_pertussis_path <- "month_sex_pertussis_20260724.csv"

# Location where you want to export the 5-year annual incidence figure
incidence_fig_path <- "fig11-2.png"


#### Setup ####

library(this.path)
library(tidyverse)
# library(slider)  # This doesn't work with every version of R, apparently.

# Import MDSS query
mdss <- read_csv(here(mdss_path), skip = 1) |> 
  
  # Clean up variable names
  janitor::clean_names() |> 
  
  # Reformat all date variables
  mutate(across(matches("date"), mdy))

# Import Ingham County population
pop <- read_csv(here(population_path)) |> select(-geo)


#### 5-year pertussis incidence per year by sex ####

# Create a table of all unique year and sex values
year_sex <- expand_grid(
  year = year(min(mdss$diagnosis_date)):year(max(mdss$diagnosis_date)),
  sex = unique(mdss$sex)
)

# Generating a 5-year sliding annual incidence rate
sliding_5yr_incidence_rate <- mdss |> 
  
  # Extract just the year of each diagnosis date
  mutate(diagnosis_year = year(diagnosis_date)) |> 
  
  # Get the year- & sex-wise cases
  summarize(cases = n(), .by = c(diagnosis_year, sex)) |> 
  
  # Use the cross-join table we made above to have every year-sex combo
  right_join(year_sex, join_by(diagnosis_year == year, sex)) |> 
  
  # Replace the missing values in the 'cases' column with 0s
  mutate(cases = cases |> replace_na(0)) |> 
  
  # Sort it, just to make it look nice
  arrange(diagnosis_year, sex) |> 
  
  # Join the populations by year & sex
  full_join(pop, join_by(diagnosis_year == year, sex)) |>
  
  # Generate 5-year rolling sum
  # Some participants may have older versions of R that don't work with slider (this happened to me hehe! - M)
  group_by(sex) |> 
  mutate(
    sliding_sum = zoo::rollsum(
      cases, 5,
      align = "right",
      na.pad = TRUE
    )
  ) |>
  
  # Generate a character vector of year ranges and the annual incidence rates
  mutate(
    period = str_c(diagnosis_year - 4, "-", diagnosis_year),
    annual_rate_per_100k = (1e5 * sliding_sum / population) / 5
  ) |>
  
  # Trim off the extra years before 2019 that we no longer need
  filter(diagnosis_year >= 2019)


# Making our pretty graph, one that *actually* looks like Danielle's
fig <- sliding_5yr_incidence_rate |> 
  
  # Implement data suppression for 5-year periods with <20 cases
  mutate(
    labels = if_else(
      condition = sliding_sum < 20,
      true = "**",
      false = sprintf("%.1f", round(annual_rate_per_100k, 1))
    ),
    # alternative to replace_when() since my tidyverse is too old!
    annual_rate_per_100k = case_when(
      sliding_sum < 20 ~ 0,
      .default = annual_rate_per_100k)
    # annual_rate_per_100k = annual_rate_per_100k |> replace_when(
    #   sliding_sum < 20 ~ 0
    # )
  ) |> 
  
  # Begin ggplot expression
  ggplot(aes(x = period, y = annual_rate_per_100k, fill = sex)) +
  geom_col(position = position_dodge(0.95)) +
  
  # Use the color scale we want
  scale_fill_manual(values = c(
    "Female" = "#73a950",
    "Male" = "#00629B"
  )) +
  
  # Remove y-axis padding below the bars (AI helped me here)
  scale_y_continuous(expand = expansion(mult = c(0, 0.05))) +
  
  # Add data labels
  geom_text(
    aes(label = labels, y = annual_rate_per_100k + 0.3),
    position = position_dodge(0.95), size = 4
  ) +
  
  # Label our axes, etc
  labs(
    title = "Pertussis\n5-Year Incidence Rate by Sex and Year 2015-2024",
    subtitle = "MADE WITH SYNTHETIC DATA",
    x = "Year", y = "Annual Rate per 100,000 Population", fill = "Legend",
    caption = str_c(
      "**Incidence rates calculated from less than 20 cases are considered",
      " statistically unreliable."
    )
  ) +
  
  # Fix some aesthetic features
  theme_classic() +
  theme(
    plot.title = element_text(hjust = 0.5),
    plot.subtitle = element_text(hjust = 0.5),
    plot.caption = element_text(hjust = 0),
    legend.position = "bottom"
  )

fig

# Export our figure
ggsave(incidence_fig_path, fig, width = 7, height = 5)


#### Cases per month by sex ####

# Generate a table of all months (as in month-year combos) that occur over the
#   period covered by 'mdss'

# Take advantage of the fact that dates are stored as integers
mdss$diagnosis_date |> head()
mdss$diagnosis_date |> as.integer() |> head()

# Get the minimum and maximum dates of our data
min_date <- min(mdss$diagnosis_date)
max_date <- max(mdss$diagnosis_date)

# Create a table with one column: all the integers including and between them
all_months <- tibble(
  # Make sure to convert the integers back to dates
  month = as_date(min_date:max_date)
) |> 
  
  # Only keep rows that represent the first day of a month
  filter(mday(month) == 1)

# Create a cross-join of all unique month and sex values
month_sex <- cross_join(
  all_months,
  tibble(sex = unique(mdss$sex))
)

# Generate a month-sex pertussis cases table
pertussis_month_sex <- mdss |> 
  mutate(diagnosis_month = diagnosis_date |> floor_date("month")) |> 
  summarize(cases = n(), .by = c(diagnosis_month, sex)) |> 
  right_join(month_sex, join_by(diagnosis_month == month, sex)) |> 
  arrange(diagnosis_month, sex) |> 
  pivot_wider(names_from = sex, values_from = cases) |> 
  mutate(across(Female:Male, \(col) col |> replace_na(0)))

# Export this table to a CSV
pertussis_month_sex |> 
  mutate(diagnosis_month = str_c(
    month(diagnosis_month), "/", year(diagnosis_month)
  )) |> 
  rename(`Diagnosis Month` = diagnosis_month) |> 
  write_csv(here(month_sex_pertussis_path))