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

# This program file demonstrates strategies discussed in
# session 4 of the 2026 NDACAN Summer Training Series 
# "Handling missing data." 

# 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. INVESTIGATING MISSING DATA
# 2. HANDLING MISSING DATA

############
# 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, 
               ggstance,
               mice, ggmice) # new missing data packages

# 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.
# Especially important with mice package. 
set.seed(1013) 

## READING DATA ## 

# Let's read in our cleaned 
# anonymized versions of the
# NCANDS Child Files for 2020-2024 (see session 2).
nc <- read_rds(paste0(data,'ncands_clean.rds'))

# Let's also 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'))

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

#################################
# 1. INVESTIGATING MISSING DATA #
#################################

## RECORD-LEVEL MISSINGNESS ## 

# Identifying record-level missingness can take elbow grease. 
# As a starting point, it's usually helpful to plot 
# record counts by reporting unit (usually state) 
# and reporting period (usually submission year). 
nc |> 
  count(staterr, subyr) |> 
  ggplot(aes(x = subyr, y = n)) + 
  geom_point() + 
  geom_line() + 
  facet_wrap(~staterr, scales = 'free') + 
  theme(axis.text.x = element_text(angle = 90))

# Recall that using Child Files as a proxy for 
# fiscal years will yield missing records toward the 
# end of the FY due to delayed reporting. 
nc |> 
  filter(rptdt %between% 
         c('2019-10-01','2020-9-30')) |> 
  group_by(rptdt, subyr) %>%
  summarise(n = n(), .groups = 'keep') |> 
  ungroup() |> 
  ggplot(aes(x = rptdt, y = n, fill = fct_rev(as.factor(subyr)))) + 
  geom_col(position = 'stack') +
  labs(fill = 'Submission year', x = 'Date', y = 'Sampled records') +
  scale_x_date(date_breaks = "1 month", date_labels = "%b %Y") +
  theme(axis.text.x = element_text(angle = 90, 
                                   hjust = .5, vjust = .5))

## ITEM-LEVEL MISSINGNESS ## 

# Item-level missingness is a little more straightforward.
# Let's quickly search for missing values 
# in our NCANDS Child Files:
# What percentage of each variable has missing values?
nc |> 
  summarize_all(~sum(is.na(.))) |> 
  mutate_all(~round(./nrow(nc)*100,2)) |> 
  t()

# Recall that searching for NA values only works
# if you've already properly cleaned your data (see session 2). 
# Don't trust the codebook blindly: 
# verify all variable encodings yourself. 
nc |> 
  count(per1rel) # original: missing data incorrectly encoded
nc |> 
  count(per1rel_new) # clean encoding

## EXPLORING THE MISSING-DATA MECHANISM ##

# As discussed, the missing-data mechanism matters. 
# While it's not usually possible to observe it directly
# from the data alone, we can make inferences about it
# using missing-data patterns. 

# For example, plotting variable distributions by 
# state and year can help identify how much missingness 
# is a function of the state or state-year of report. 
nc |> 
  count(staterr, subyr, fcmoney) |> 
  group_by(staterr, subyr) |> 
  mutate(p = n/sum(n)) |> 
  ggplot(aes(x = p, y = staterr, fill = fcmoney)) + 
  geom_histogram(stat = 'identity') + 
  facet_wrap(~ subyr) +
  theme(axis.text.x = element_text(angle = 90, 
                                   vjust = .5))

# This gives us good reason to suspect that most 
# missingness for the fcmoney variable results from 
# (changes in) state-level data collection and 
# reporting protocols. It also suggests that values 
# for fcmoney are (mostly) missing at random rather than
# missing completely at random or missing not at random. 

# Note that in the case of linked data, 
# *item*-level missingnesscan arise from 
# *record*-level missingness.
# Examining our linked dataset, we see that a 
# (very small) number of observations have missing values
# for maltreatment history because these children had 
# missing records of maltreatment.
dlink |> 
  summarize_all(~sum(is.na(.))) |> 
  mutate_all(~round(./nrow(nc)*100,2)) |> 
  t()

############################
# 2. HANDLING MISSING DATA #
############################

# Let's explore options for handling missing data 
# using maltreatment reports for FY2020-2024. 
# Let's also create an indicator variable measuring
# whether a report was substantiated or indicated.
ncs <- nc |> 
  filter(rptdt %between% 
           c('2019-10-01','2023-9-30')) |>
  mutate(subind = case_when(
           rptdisp %in% c('Substantiated',
                          'Indicated/reason to suspect') ~ 1,
           is.na(rptdisp) ~ NA, 
           T ~ 0)
         )

# COMPLETE-CASE ANALYSIS # 

# Let's estimate a very basic logistic regression model, 
# and summarize it. Counting the observations used in the model, 
# note that it differs from the length of the full dataset. 
# This is because most models can't be estimated directly on 
# observations with missing values of modeled variables. 
# In other words, the default is to conduct a 
# complete-case analysis, or to listwise-delete 
# observations with missing values.
m_cc <- glm(subind ~ chprior + fcmoney, 
            data = ncs, 
            family = 'binomial') 
summary(m_cc)
nobs(m_cc)
nrow(ncs)

# But listwise deletion is defensible ONLY if our data are MCAR. 
# Recall that our data are most likely mostly MAR.

# HOT DECK IMPUTATION #

# Hot deck imputations fills in missing values 
# using observed values from other units with 
# similar observed values. Usually it's not a very good way 
# to impute missing data. But the NCANDS Child Files 
# (and other NDACAN data) have a special property that makes it 
# suitable: repeated observations of the same units, in this case, 
# multiple reports for many of the same children. Perhaps fcmoney
# is missing for a child on one report but not on another. 
# So let's sort by child ID and then report date.
ncs |> 
  arrange(stfcid, rptdt) 

# And count rows representing the same child and for which 
# fcmoney is missing in one row but not another.
ncs |> 
  filter(is.na(fcmoney) & !is.na(lag(fcmoney)) & 
           stfcid == lag(stfcid) & 
           chid != 'XXXXXXXXXXXX')

# Bingo! Just one, but not nothing, and likely to be more common 
# depending on the sample and variables analyzed. 
# So let's carry over the non-missing values
# to the missing ones. Validity rests on the assumption that
# children's true financial distress is time-invariant. 
# (This assumption is actually false, 
# but might not be in other cases.) 
ncs <- ncs |> 
  group_by(stfcid) |> 
  fill(fcmoney, .direction = 'downup') |> 
  ungroup()

# MULTIPLE IMPUTATION #

# Multiple imputation is a model-based strategy that uses 
# information about the relationships between observed values
# to make guesses about the true value of missing values. 
# We make multiple guesses, or imputations, the variance of which
# captures our uncertainty about the accuracy of our model. 

# For demonstration purposes, let's take a random 3% sample
# or our data, keeping only those variables of interest. 
# Note that multiple imputation is computationally intensive. 
# Plan to have your machine impute your data while you 
# eat lunch or sleep. Consider the need for non-local computing.  
ncsi <- ncs |> 
  select(subyr, staterr, rptdt, stfcid,
         subind, 
         chprior, fcmoney) |> 
  mutate(chprior = factor(chprior, 
                          levels = c('No', 'Yes'), 
                          labels = c('No', 'Yes')), 
         fcmoney = factor(fcmoney, 
                          levels = c('No', 'Yes'), 
                          labels = c('No', 'Yes'))) |> 
  slice_sample(prop = .03)

# We examine patterns of missingness in the sample, 
# which is important for assessing the feasibility 
# of multiple imputation. 
ncsi |> 
  plot_pattern(rotate = T) + 
  theme(axis.text.x = element_text(vjust = 0.25), 
        legend.position = 'none') 

# And designate variables we *don't* want to use in the imputation.
pred <- quickpred(ncsi, 
                  exclude = c("stfcid", "rptdt"))

# Then we estimate an imputation model. 
imp <- mice(ncsi, 
            predictorMatrix = pred, 
            m = 5, 
            maxit = 5, 
            print = F)

# We can check basic features of the imputation. 
# Note that we probably wanted mice to use "logreg" for subind
# instead of "pmm." We can either recode the variable, 
# or specify imputation methods using the "method" option of mice.
summary(imp)

# We can easily diagnose the convergence of our imputation model 
# (see additional resources).  
plot(imp)

# There is *much* more to learn about mice. 
# For example, we may want to predict missing values
# based not just on state and year, but based on state-years. 
# At very least explore the help documentation, 
# including helpful vignettes. 
help(mice)

# Let's now use our imputed data to estimate the same, 
# basic logistic regression model that we estimated 
# on our complete-case data. 
m_mice <- with(imp, 
               glm(subind ~ chprior + fcmoney,
                   family = 'binomial'))

# And visually compare the results for the two models. 
# First, re-estimate our complete-count model on the same
# random 3% sample. 
m_cc2 <- glm(subind ~ chprior + fcmoney, 
            data = ncsi, 
            family = 'binomial') 

# Then organize the estimates. 
cc_summary <- summary(m_cc2) |>
  coef() |>
  as.data.frame() |>
  rownames_to_column('term') |>
  mutate(model = 'CC') |>
  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')

# Organize the mice estimates. 
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)

# And finally, combine and visually compare the estimates. 
bind_rows(mice_combined, cc_combined) |>
  mutate(term = factor(term, 
                       levels = c('(Intercept)', 
                                  'chpriorYes', 
                                  'fcmoneyYes')), 
         est   = exp(est),
         lower = exp(lower),
         upper = exp(upper)) |>
  filter(term != '(Intercept)') |> 
  ggplot(aes(x = est, y = fct_rev(term),
             xmin = lower, xmax = upper,
             color = model, group = model)) +
  geom_vline(xintercept = 1, linetype = 'dashed') +
  geom_point(position = position_dodgev(height = -.5)) +
  geom_errorbarh(height = .25,
                 position = position_dodgev(height = -.5)) +
  labs(x = 'Odds ratio', y = 'Variable', color = 'Model')

# Notice that the differences are meaningful. A reminder that
# these data are fabricated, but illustrate the main point: 
# improper handling of missing data can meaningfully bias 
# your analysis and lead you to false inferences.