## Quarto

Quarto enables you to weave together content and executable code into a finished presentation. To learn more about Quarto presentations see <https://quarto.org/docs/presentations/>.

## Bullets

When you click the **Render** button a document will be generated that includes:

-   Content authored with markdown
-   Output from executable code

## Code

When you click the **Render** button a presentation will be generated that includes both content and the output of embedded code. You can embed code like this:

```{r}
#! echo: false

# 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) # to be consistent with how they'll do it (also works better on my machine anyhow for whatever reason!)
# here::i_am("pertussis_analysis.R")
# 
# library(here)
library(tidyverse)
#library(slider)

# 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 cross-join of all unique year and sex values
year_sex <- cross_join(
  # Mimicking the cross-join in the previous section in fewer lines
  tibble(year = year(min(mdss$diagnosis_date)):year(max(mdss$diagnosis_date))),
  tibble(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) |>
  
  # 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))
    ),
    annual_rate_per_100k = case_when(sliding_sum < 20 ~ 0, .default = annual_rate_per_100k) # alternative since my tidyverse is also too old for replace_when!
    # annual_rate_per_100k = annual_rate_per_100k |> replace_when(
    #  sliding_sum < 20 ~ 0
    #)
  ) 


# Making our pretty graph, one that *actually* looks like Danielle's
fig <- sliding_5yr_incidence_rate |> 
  
  # 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
```


## Lollipop

```{r}
# Library
library(ggplot2)
library(dplyr)
# library(hrbrthemes)

# Create data
value1 <- abs(rnorm(26))*2
data <- data.frame(
  x=LETTERS[1:26], 
  value1=value1, 
  value2=value1+1+rnorm(26, sd=1) 
)
 
# Reorder data using average? Learn more about reordering in chart #267
data <- data %>% 
  rowwise() %>% 
  mutate( mymean = mean(c(value1,value2) )) %>% 
  arrange(mymean) %>% 
  mutate(x=factor(x, x))
 
# Plot
ggplot(data) +
  geom_segment( aes(x=x, xend=x, y=value1, yend=value2), color="grey") +
  geom_point( aes(x=x, y=value1), color=rgb(0.2,0.7,0.1,0.5), size=3 ) +
  geom_point( aes(x=x, y=value2), color=rgb(0.7,0.2,0.1,0.5), size=3 ) +
  coord_flip()+
  theme(
    legend.position = "none",
  ) +
  xlab("") +
  ylab("Value of Y")
```