# intro_code_ref.R
#
# Script containing all code used to demonstrate R to attendees of MI-Support's
# 24 July 2026 Intro to R workshop
#
# Author: Jack Jacobs (jckjcbs@umich.edu)
#         MI-Support (micom-misupport@umich.edu)

#### Setup ####

# install.packages(c("tidyverse","this.path","here","janitor"))

# Load tidyverse functions for data analysis
library(tidyverse)

# Ask your computer "where you are" for this exercise
this.path::this.path()

# Example of how 'this.path' can help you
this.path::here("fake file name")


#### Importing data ####

# "Read in" fake data about pertussis into your session
mdss <- read_csv(here("data/FAKE_pertussis_mdss.csv"), skip = 1)
View(mdss)

# Let's clean this table up a bit.
# Our first pipeline: remove some columns and clean the names
mdss <- mdss |> 
  
  # Grab just the columns we care about
  select(
    InvestigationID, Onset_Date, Diagnosis_Date, Referral_Date, Patient_Status,
    Case_Disposition, Date_of_Death, Date_of_Birth, Sex, Age, Admission_Date
  ) |> 
  
  # janitor is a nice little package
  # https://sfirke.github.io/janitor/articles/janitor.html
  janitor::clean_names()

# Quick preview
mdss |> 
  ggplot(aes(x = year(mdy(diagnosis_date)), fill = sex)) +
  geom_col(stat = "count", position = "dodge")


#### Backing up, what *is* R? ####

# So... what is this data? What's happening here?
str(mdss)
mdss$investigation_id
nrow(mdss)
dim(mdss)
mdss$age
length(mdss$age)
unique(mdss$age)
mdss$age |> unique()
n_distinct(mdss$age)
table(mdss$age)
colnames(mdss) 

# Let's back up even further... what is R??
# (A programming language: something to use to tell your computer what to do)
cat("Hello world!")
3 + 7

# R is often called a "statistical programming language"
# Here are some basic operations
c(3 + 7, 12.6 + 1.234, 0.9876 / 2)
numbers <- c(3 + 7, 12.6 + 1.234, 0.9876 / 2)
round(numbers)
mean(numbers)
is.na(numbers)  # Useful when you have missing data

# There are other data types: here are "characters" or "strings"
c("hello", "hi", "hey there!")
c("hello") == "hello"
"hello" == "hi"

# And here are "logicals"
TRUE & FALSE
TRUE | FALSE
!TRUE
logicals <- c(T, T, F)
any(logicals)
all(logicals)

# Strings can be kindof strange
"hello" < "hi"
spaced <- "1 2 3 4 5 6 7 8 9 10"
spaced |> str_split_1(" ") |> as.numeric()
A_spaced <- paste0("A", str_split_1(spaced, " "))
print(A_spaced)
print(A_spaced[3])
print(A_spaced[2:6])
2:6

# But their flexibility can be a good thing!
R_birthday <- "8/16/1993"
R_birthday <- R_birthday |> mdy()
class(R_birthday)
R_age <- R_birthday %--% today()
as.period(R_age)


#### Data cleaning ####

# Read in a population table. We'll use this in a minute
pop <- read_csv(here("population.csv"))

# Making more changes to our data in a pipeline with lots of comments!
mdss <- mdss |>
  
  # Making some useful variables with mutate
  mutate(
    # Use 'lubridate' to make a date vector out of date-like strings
    diagnosis_date = mdy(diagnosis_date),
    
    # You can immediately reference variables you make within a mutate function
    diagnosis_year = year(diagnosis_date)
  ) |> 
  
  # Use filter to remove rows based on some T/F condition(s)
  filter(diagnosis_year >= 2019) |> 
  
  # Renaming is easy
  rename(
    # new_name = old_name
    age_years = age
  )

# Let's find some aggregate variables
summary <- mdss |> 
  group_by(diagnosis_year, sex) |> 
  summarize(cases = n(), .groups = "drop")


##### Joins #####

# Trying to put two tables together by lining up values? That's a join.

# There weren't any cases among Females in 2020!
# Ensure you have all categories by "expanding the grid" of your options
all_categories <- expand_grid(
  diagnosis_year = 2019:2024,
  sex = c("Female","Male")
)

# Join these categories to your data (show Venn diagram for "right-join")
summary <- right_join(
  summary, all_categories,
  join_by(diagnosis_year, sex)
) |> 
  replace_na(list(cases = 0)) |>  # Replace the missing cases value with 0
  arrange(diagnosis_year, sex)

# We want to incorporate the scale of each state's population!
# (show Venn diagram for "left-join")
summary <- left_join(
  summary, pop,
  # by = c("diagnosis_year" = "year", "sex")
  join_by(diagnosis_year == year, sex)
)

# Now that we have the population, we can calculate rates! (per 100k pop)
viz_df <- summary |> 
  mutate(cases_per_100k = 1e5 * cases / population)


#### Making graphs ####

# Before we get going, we should acknowledge that this is... a lot. It's a lot
#   of detail and fine-tuning, and it takes a while to get right. However,
#   **it's one and done!** Excel is similar in the level of messing with 
#   settings required, but this way, you can "transport" this to similar graphs
#   you make in the future. In Excel, you need to start mostly from scratch for
#   each new graph.

# Let's look at what our bar graph will roughly look like
viz_df |> 
  ggplot(aes(x = diagnosis_year, y = cases_per_100k, fill = sex)) +
  geom_col(position = "dodge")

# Let's align the aesthetics with our reference graph
viz_df |> 
  ggplot(aes(x = diagnosis_year, y = cases_per_100k, fill = sex)) +
  geom_col(position = "dodge") +
  scale_fill_manual(
    "Sex",
    values = c("Female" = "#73A950", "Male" = "#00629B")
  ) +
  theme_classic()

# Axis labels & title
viz_df |> 
  mutate(diagnosis_year = as.character(diagnosis_year)) |> 
  ggplot(aes(x = diagnosis_year, y = cases_per_100k, fill = sex)) +
  geom_col(position = "dodge") +
  scale_fill_manual(
    "Sex",
    values = c("Female" = "#73A950", "Male" = "#00629B")
  ) +
  labs(
    title = "Pertussis\nIncidence Rate by Sex and Year 2019-2024",
    subtitle = "Ingham County",
    x = "Year", y = "Rate per 100,000 Population"
  ) +
  theme_classic() +
  theme(
    plot.title = element_text(hjust = 0.5),
    plot.subtitle = element_text(hjust = 0.5),
    legend.position = "bottom"
  )

# Let's add quantity labels!
viz_df |> 
  mutate(
    diagnosis_year = as.character(diagnosis_year),
    text = cases_per_100k |> round(1) |> as.character()
  ) |> 
  ggplot(aes(x = diagnosis_year, y = cases_per_100k, fill = sex)) +
  geom_col(position = "dodge") +
  scale_fill_manual(
    "Sex",
    values = c("Female" = "#73A950", "Male" = "#00629B")
  ) +
  geom_text(
    aes(label = text, y = cases_per_100k + 1.2),
    position = position_dodge(0.9),
    size = 3
  ) +
  labs(
    title = "Pertussis\nIncidence Rate by Sex and Year 2019-2024",
    subtitle = "Ingham County",
    x = "Year", y = "Rate per 100,000 Population"
  ) +
  theme_classic() +
  theme(
    plot.title = element_text(hjust = 0.5),
    plot.subtitle = element_text(hjust = 0.5),
    legend.position = "bottom"
  )

# Data suppression (for demonstration, we'll suppress below counts of 3)
viz_df |> 
  mutate(
    diagnosis_year = as.character(diagnosis_year),
    text = cases_per_100k |> round(1) |> as.character(),
    text = if_else(
      cases < 3, true = "**", false = text
    ),
    cases_per_100k = if_else(
      cases < 3, true = 0, false = cases_per_100k
    )
  ) |> 
  ggplot(aes(x = diagnosis_year, y = cases_per_100k, fill = sex)) +
  geom_col(position = position_dodge(0.95)) +
  scale_fill_manual(
    "Sex",
    values = c("Female" = "#73A950", "Male" = "#00629B")
  ) +
  geom_text(
    aes(label = text, y = cases_per_100k + 1.2),
    position = position_dodge(0.95),
    size = 3
  ) +
  labs(
    title = "Pertussis\nIncidence Rate by Sex and Year 2019-2024",
    subtitle = "Ingham County",
    x = "Year", y = "Rate per 100,000 Population",
    caption = str_c(
      "**Incidence rates calculated from less than 3 cases are considered ",
      "unreliable."
    )
  ) +
  theme_classic() +
  theme(
    plot.title = element_text(hjust = 0.5),
    plot.subtitle = element_text(hjust = 0.5),
    legend.position = "bottom",
    plot.caption = element_text(hjust = 0)
  )

# Finally, let's save our plot
ggsave(here("img/our_final_plot.png"))