R version 4.5.2 (2025-10-31 ucrt) -- "[Not] Part in a Rumble" Copyright (C) 2025 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. Natural language support but running in an English locale R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > # 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) > # 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") > # 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() Warning message: Removed 189 rows containing missing values or values outside the scale range (`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 Warning message: Removed 189 rows containing missing values or values outside the scale range (`geom_point()`). > p1 + + geom_point(position = position_jitter(width = .2, + height = .2)) # place noisily Warning message: Removed 189 rows containing missing values or values outside the scale range (`geom_point()`). > # 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() `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")' Warning messages: 1: Removed 189 rows containing non-finite outside the scale range (`stat_smooth()`). 2: Removed 189 rows containing missing values or values outside the scale range (`geom_point()`). > # 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')) `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")' Warning messages: 1: Removed 189 rows containing non-finite outside the scale range (`stat_smooth()`). 2: Removed 189 rows containing missing values or values outside the scale range (`geom_point()`). > # 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 # A tibble: 2 × 6 var Mean Median SD P10 P90 1 Substantiated investigations 1.36 1 0.877 1 3 2 Age 0.827 0 1.05 0 3 > # 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) \begin{tabular}{lrrrrr} \toprule var & Mean & Median & SD & P10 & P90\\ \midrule Substantiated investigations & 1.3648000 & 1 & 0.8772849 & 1 & 3\\ Age & 0.8267014 & 0 & 1.0492420 & 0 & 3\\ \bottomrule \end{tabular} > # 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') Warning message: Removed 189 rows containing non-finite outside the scale range (`stat_bin()`). > # 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 strenghten 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 %`) Waiting for profiling to be done... > 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') Warning messages: 1: position_dodgev requires non-overlapping y intervals 2: Using the `size` aesthetic with geom_path was deprecated in ggplot2 3.4.0. ℹ Please use the `linewidth` aesthetic instead. This warning is displayed once per session. Call lifecycle::last_lifecycle_warnings() to see where this warning was generated.