#########
# NOTES #
#########

# This program file demonstrates strategies discussed in
# session 5 of the 2026 NDACAN Summer Training Series 
# "Data presentation and visualization." 

# For questions, contact the presenter
# Alex Roehrkasse (aroehrkasse@butler.edu).

# Note that because of the process used to anonymize data, 
# all unique observations include partially fabricated data
# that prevent the identification of respondents.
# As a result, all descriptive and 
# model-based results are fabricated.

# Results from this and all NDACAN presentations are for 
# training purposes only and should never be understood 
# or cited as analysis of NDACAN data.

#####################
# TABLE OF CONTENTS #
#####################

# 0. SETUP
# 1. INTRO TO GGPLOT2
# 2. DESCRIPTIVE STATISTICS
# 3. MODEL STATISTICS

############
# 0. SETUP #
############

## SETTING UP THE ENVIRONMENT ## 

# Let's clear the environment.
rm(list=ls())

# Pacman installs packages if necessary, otherwise loading them.
if (!requireNamespace("pacman", quietly = TRUE)){
  install.packages("pacman")
}
pacman::p_load(data.table, tidyverse, 
               mice,
               knitr,
               ggstance, gtsummary, modelsummary, scales)

# Let's define some filepaths (note the organization of project and data folders).
project <- 'C:/Users/aroehrkasse/Box/Presentations/-NDACAN/2026_summer_series/'
data <- 'C:/Users/aroehrkasse/Box/NDACAN/2026_summer_series/'

# And set one as the working directory.
setwd(project)

# Always set a seed to allow for reproduction of random processes.
set.seed(1013) 

## READING DATA ## 

# Let's read in our cleaned, linked data: children 0-3 entering foster care in 2023 (AFCARS) linked to maltreatment histories (NCANDS) (see session 3). 
dlink <- read_rds(paste0(data,'linked_data.rds'))

# Let's also load our complete-case and multiply imputed model estimates of substantiation/indication (see session 4). 
load("saved_models.RData")

# For the purposes of the presentation, note again that I have sampled New England and selected a small number of key variables. 


#######################
# 1. INTRO TO GGPLOT2 #
#######################

# The modern grammar of data visualization in R is implemented in the ggplot2 package, part of the Tidyverse. 
# To create a ggplot, you first apply the function to a dataset and define aesthetics. Notice, though, that we still don't have any data visualized. 
p1 <- ggplot(dlink, 
             aes(x = nrep,
                 y = nsub)) 
p1 

# This is because we then have to add geoms. 
p1 + 
  geom_point()

# We can then modify these geoms by adding arguments to them. Note how each approach (imperfectly) leverages pre-attentive attributes. 
p1 + 
  geom_point(alpha = .05) # set transparency
p1 + 
  geom_point(position = position_jitter(width = .2, 
                                        height = .2)) # place noisily

# Furthermore, we can overlay multiple geoms. 
p1 + 
  geom_point(position = position_jitter(width = .2, 
                                        height = .2)) + 
  # A line illustrating identity between the two axes
  geom_abline(slope = 1, 
              linewidth = 1,
              linetype = 'dashed') + # accessible encoding
  # A locally estimated line of best fit
  geom_smooth()

# But to really know what's going on, we need to label things clearly and completely. 
p1 + 
  geom_point(position = position_jitter(width = .2, 
                                        height = .2)) + 
  # Make color an aesthetic so it can be included in a legend
  geom_abline(aes(intercept = 0, 
                  slope = 1, 
                  color = '100% substantiation'), 
              linewidth = 1, 
              linetype = 'dashed') + 
  geom_smooth(aes(color = 'Locally fitted line')) + 
  # Set colors manually
  scale_color_manual(values = c('100% substantiation' = 'red', 
                                'Locally fitted line' = 'blue')) + 
  # Label axes and suppress legend title
  labs(x = 'Number of maltreatment reports', 
       y = 'Number of\nsubstantiated reports', 
       color = NULL) + 
  # Choose theme
  theme_classic() + 
  # Design legend
  theme(legend.position = c(0.01, 0.99),
        legend.justification = c("left", "top"), 
        legend.background = element_rect(color = 'black'))

# The figure, while far from perfect, communicates a few things clearly and intuitively: (1) most children had a small number of reports and substantiated reports, but (2) some had many, and (3) the proportion of reports that were substantiated decreased with the number of reports.   

#############################
# 2. DESCRIPTIVE STATISTICS #
#############################

# Very conventional in an academic journal article is a table of descriptive statistics. Take, for example, our linked dataset. Some helpful canned functions can generate nice-looking descriptive tables, e.g. the gtsummary package.
dlink |> 
  tbl_summary(include = c(nsub, ageatlatrem))

# Sometimes, though, manual calculation is more flexible.
dlink_sum <- dlink |> 
  select(nsub, ageatlatrem) |> 
  summarize(across(everything(), 
                   list(
                     Mean = ~mean(., na.rm = T),
                     Median = ~median(., na.rm = T),
                     SD = ~sd(., na.rm = TRUE),
                     P10 = ~quantile(., 0.1, na.rm = T),
                     P90 = ~quantile(., 0.9, na.rm = T)
                   ))) |> 
  pivot_longer(cols = everything(), 
               names_to = c('var', '.value'), 
               names_sep = '_') |> 
  mutate(var = case_when(var == 'ageatlatrem' ~ 'Age', 
                         var == 'nsub' ~ 'Substantiated investigations'))
dlink_sum

# We can save this table as a CSV or XLSX file, or output it as a table in LaTeX using the knitr package.  
kable(dlink_sum, format = "latex", booktabs = TRUE)

# Consider, however, how visualization can combine some of the merits of each of these types of tables, intuitively presenting both the full distribution of the variables and their central tendencies.  
dlink |> 
  select(nsub, ageatlatrem) |> 
  pivot_longer(cols = everything(), 
               names_to = 'var', values_to = 'num') |> 
  mutate(var = case_when(var == 'ageatlatrem' ~ 'Age', 
                         var == 'nsub' ~ 'Substantiated investigations')) |>
  ggplot() +
  # Plot a histogram with width of 1
  geom_histogram(aes(x = num),
                 binwidth = 1) + 
  # Use our summarized data from above to plot central tendencies
  geom_point(data = dlink_sum |> 
               pivot_longer(cols = c('Mean', 'Median'), 
                            names_to = 'stat'), 
             aes(x = value, 
                 # Multiply encode the central tendencies
                 color = stat, 
                 shape = stat), 
             y = 0, 
             size = 3) +
  # Design scales
  scale_x_continuous(breaks = 0:10) + 
  scale_y_continuous(labels = label_comma()) + 
  # Choose an accessible color palette
  scale_color_brewer(palette = 'Set2') + 
  # Create separate panels for each variable
  facet_wrap(~var, scales = 'free_x') + 
  theme_bw() +
  guides(color = guide_legend(reverse = T), 
         shape = guide_legend(reverse = T)) + 
  labs(x = NULL, y = 'Number of\nobservations', 
       color = NULL, shape = NULL) + 
  theme(legend.position = 'bottom')

#######################
# 3. MODEL STATISTICS #
#######################

# Recall from session 3 that we estimated a basic logistic regression model using complete-case analysis and multiple imputation. Conventionally, we would report these results as a table. Again, there are some helpful canned functions for table generation (cf. the stargazer package). 
modelsummary(list(m_cc2, pool(m_mice)))

# These packages make it easy to customize your tables to add/remove/label information to make your message clearer. 
modelsummary(list('Complete case' = m_cc2, 
                  'MICE' = pool(m_mice)),
             statistic = "conf.int", 
             conf_level = 0.95, 
             stars = T, 
             coef_omit = "Intercept", 
             coef_rename = c("chpriorYes" = "Prior report", 
                             "fcmoneyYes" = "Financial distress"),
             gof_omit = ".*")

# While this table is clear, it's not altogether intuitive. Consider how visualization can strengthen our evidence-based message. 
# First, we extract, clean, and combine quantities of interest from our models, creating a plot_data data frame.
cc_summary <- summary(m_cc2) |>
  coef() |>
  as.data.frame() |>
  rownames_to_column('term') |>
  mutate(model = 'Complete case') |>
  rename(est = Estimate) |>
  select(term, est, model)
cc_ci <- confint(m_cc2) |>
  as.data.frame() |>
  rownames_to_column('term') |>
  rename(lower = `2.5 %`, upper = `97.5 %`)
cc_combined <- left_join(cc_summary, cc_ci, by = 'term')
mice_combined <- pool(m_mice) |>
  summary(conf.int = TRUE) |>
  as.data.frame() |>
  mutate(model = 'MICE') |>
  rename(est = estimate, lower = `2.5 %`, upper = `97.5 %`) |>
  select(term, est, lower, upper, model)
plot_data <- bind_rows(mice_combined, cc_combined) |>
  mutate(term = factor(term, 
                       levels = c('(Intercept)', 
                                  'chpriorYes', 
                                  'fcmoneyYes'), 
                       labels = c('Intercept', 
                                  'Prior report', 
                                  'Financial distress')), 
         est   = exp(est),
         lower = exp(lower),
         upper = exp(upper)) |>
  filter(term != 'Intercept') 

# We can then plot the quantities in the plot_data data frame.  
plot_data |> 
  # Define aesthetics shared across geoms
  ggplot(aes(x = est, 
             y = fct_rev(term),
             color = model, 
             # Multiply encode for accessibility
             shape = model, 
             group = model)) +
  # Highlight the null hypothesis as a visual anchor
  geom_vline(xintercept = 1, linetype = 'dashed') + 
  # Plot point estimates as points
  geom_point(position = position_dodgev(height = -.5), 
             size = 2) + 
  # Plot confidence intervals as error bars
  geom_errorbarh(aes(xmin = lower, 
                     xmax = upper), 
                 height = .25,
                 position = position_dodgev(height = -.5)) +
  # Log-transform axis of ratio estimand for comparison above/below 1
  scale_x_continuous(trans = 'log', 
                     breaks = seq(.5,2.25,.25), 
                     limits = c(.5, 2.25)) + 
  # Choose a colorblind-safe color palette 
  scale_color_brewer(palette = 'Dark2') + 
  # Clearly label all plot features
  labs(x = 'Odds ratio', y = 'Predictor', 
       color = 'Missing data\nstrategy', 
       shape = 'Missing data\nstrategy') + 
  ggtitle('Substantiation/indication of\nCPS reports') + 
  theme_bw() + 
  theme(axis.text.x = element_text(angle = 45, 
                                   hjust = 1, 
                                   vjust = 1.1), 
        legend.position = 'bottom')