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

# This program file demonstrates strategies discussed in
# session 3 of the 2026 NDACAN Summer Training Series 
# "Linking NCANDS and AFCARS." 

# 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. LINKING BY STACKING: AGGREGATE DATA
# 2. LINKING BY STACKING: MICRODATA
# 3. LINKING BY JOINING: MICRODATA


############
# 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) 

# 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 
# anonymized versions of the
# NCANDS Child Files for 2020-2024.
nc <- read_rds(paste0(data,'ncands_clean.rds'))

# Let's also read in cleaned 
# anonymized versions of the 
# AFCARS Foster Care AB files for 2023 and 2024.
ac23 <- read_rds(paste0(data,'afcars23_clean.rds'))
ac24 <- read_rds(paste0(data,'afcars24_clean.rds'))

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

##########################################
# 1. LINKING BY STACKING: AGGREGATE DATA #
##########################################

# Recall from S2 that nc is already "stacked" 
# after we unlisted our NCANDS submission-year 
# into a single data frame. 
nc |> 
  count(staterr,subyr)

# We can do the same explicitly by row-binding
# multiple AFCARS AB files. 
# Note that this requires consistent variable naming and encoding.
ac <- ac23 |> 
  bind_rows(ac24)

ac |> 
  count(st, fy)

# Say we wanted to link states across years, 
# measuring the annual proportion of maltreatment reports 
# that were substantiated. 
# This requires using NCANDS to link at the state-year level 
# by stacking and summarizing. 
# Recall also that we need Y+1 years of data to avoid
# bias from delayed reporting, so let's examine reports
# from FY2020-2023 using data from FY2020-2024.
state_year <- nc |> 
  filter(rptdt %between% 
           c('2019-10-01','2023-9-30')) |> # define sampling frame
  mutate(fy = if_else(month(rptdt) >= 10, 
                      year(rptdt) + 1, 
                      year(rptdt)), # create fiscal year variable
         sub = case_when(rptdisp == 'Substantiated' ~ 1, 
                         is.na(rptdisp) ~ NA, 
                         T ~ 0)) |> # create substantiation indicator
  group_by(fy, staterr) |> # group by year and state
  summarize(nsub = sum(sub, na.rm = T), # count substantiated reports
            n = n(), # count all reports
            .groups = 'drop') |> 
  mutate(sub_prop = nsub/n) |> # create substantiation proportion variable
  arrange(staterr, fy) # organize

# Notice that we now have a panel of state-years
# with the linking variables (a combination of state and year)
# and the outcome of interest over time (substantiation proportion). 
state_year |> 
  ggplot(aes(x = fy, y = sub_prop, 
             color = staterr, shape = staterr)) + 
  geom_point() + 
  geom_line()


######################################
# 2. LINKING BY STACKING: MICRO DATA #
######################################

# Linking microdata by stacking requires less manipulation
# but much more care and caution. 
# Let's illustrate with the AFCARS AB files.

# Let's try to link individual children across the 
# 2023 and 2024 fiscal years. In principle, 
# once our AB files are stacked, this only requires 
# use of the linking variable, the unique child ID StFCID, 
# to identify children over time. 
# Notice that with just a little rearranging, we identify
# linked children: rows with matching values of the 
# linking variable StFCID. 
ac |> 
  arrange(stfcid, fy) |> # order descending by state, then year
  select(fy, st, stfcid, entered, exited, inatend, inatstart) |> 
  slice(1:10) # pick out the first 10 rows

# We can examine our raw linkage rates by state and year: 
# What proportion of children with records in one year
# also have a record in the other year? 
link_success <- function(df) {
  df |> 
    arrange(stfcid, fy) |> # important to do whenever you use lead()/lag()
    mutate(linked = if_else(fy == 2023 & stfcid == lead(stfcid) | 
                              fy == 2024 & stfcid == lag(stfcid), 
                            1, 
                            0)) |> # create a linkage indicator
    group_by(st, fy, linked) |> # group to count
    summarize(n = n()) |> # count
    group_by(st, fy) |> # regroup
    mutate(linked_prop = n/sum(n)) |> # calculate linkage proportion
    filter(linked == 1) |> # keep only success rates
    ggplot(aes(x = linked_prop, y = st, 
               color = factor(fy), shape = factor(fy))) + # visualize
    geom_point(position = ggstance::position_dodgev(height = .4)) + 
    scale_x_continuous(limits = c(0,1))
}

ac |> 
  link_success()

# But this doesn't tell us much, because children can fail
# to have repeat records for two completely different reasons: 
# (1) they were actually only in foster care for 1 of 2 years, or 
# (2) they were in foster care both years but we failed to link them. 
# We can use additional variables to evaluate our link quality, 
# e.g. information about whether a child was in care at the 
# beginning or end of the reporting period.
# If we limit our data to children whose records 
# we should be able to link, we can get a better sense of our 
# linkage quality.  
ac |> 
  filter((fy == 2023 & inatend == 1) | 
           (fy == 2024 & inatstart == 1)) |> 
  link_success()

# Another concern, however, is false positives, 
# i.e. links that don't exist but which we mistakenly observe. 
# We can use information like date of birth to verify 
# that linked children have consistent time-invariant attributes.
ac |> 
  arrange(stfcid, fy) |> 
  mutate(linked = if_else(fy == 2023 & stfcid == lead(stfcid) | 
                            fy == 2024 & stfcid == lag(stfcid), 
                          1, 
                          0), # create a linkage indicator
         dob_match = if_else(fy == 2023 & dob == lead(dob) | 
                               fy == 2024 & dob == lag(dob), 
                             1, 
                             0)) |> # create birthdate match indicator
  filter(linked == 1) |> # keep only successful links
  group_by(st, fy, dob_match) |> 
  summarize(n = n()) |> 
  group_by(st, fy) |> 
  mutate(false_pos_prop = n/sum(n)) |> # calc. prop. with matched DOB
  filter(dob_match == 1) |> # keep only success rates
  ggplot(aes(x = dob_match, y = st, 
             color = factor(fy), shape = factor(fy))) + # visualize
  geom_point(position = ggstance::position_dodgev(height = .4)) + 
  scale_x_continuous(limits = c(.9,1))

# Looks good in this toy example, but don't count on it with real data! 

######################################
# 3. LINKING BY JOINING: MICRO DATA #
######################################

# Research question: 
# What are the maltreatment histories of children 
# aged 0-3 (NCANDS) placed into foster care in FY2023 (AFCARS)?
  
## PREPARING DATA FOR ONE-TO-ONE JOIN ##

# Most essential for linking by joining one-to-one are three things: 
# 1) In each linking data set, a row corresponds to the unit of analysis. 
# 2) Any and all linking variables are identically named and encoded. 
# 3) Any and all non-linking variables are differently named. 

# First, let's use linking by stacking using NCANDS to 
# generate two maltreatment history measures: 
# total reports and total substantiated reports.
start <- Sys.time()
nc_link <- nc |> 
  filter(rptdt %between% 
           c('2019-10-01','2023-9-30')) |> # define sampling frame
  mutate(sub = case_when(rptdisp == 'Substantiated' ~ 1, 
                         is.na(rptdisp) ~ NA, 
                         T ~ 0)) |> # create substantiation indicator
  group_by(stfcid) |> # group by unique child ID
  summarize(nsub = sum(sub, na.rm = T), # count substantiated reports
            nrep = n(), # count all reports
            #maxdt = max(rptdt), # computationally intensive
            .groups = 'drop')
end <- Sys.time()
end - start
#write_rds(nc_link, paste0(data,'nc_link.rds')) 
nc_link <- read_rds(paste0(data,'nc_link.rds')) # reads pre-processed data

# Let's check the output. Something looks amiss: 
# one observation has an implausible number of reports.
# This arises from a missing value for the ID variable. 
# We'll ignore this for the purposes of demonstration, but 
# a careful analysis would inquire into the 
# causes and consequences. 
nc_link |> count(nrep) 
nc_link |> count(nsub) 

nc_link |> filter(nrep > 20 | nsub > 20)

head(nc_link)

# Let's prepare our AFCARS data. 
# We want to select only those children aged 0-3 
# entering foster care in FY2023. 
ac_link <- ac |> 
  filter(dob %between% 
           c('2019-10-01','2023-9-30') &  # define sampling frame 
           entered == 1 & 
           fy == 2023)
head(ac_link)

# Note that: 
# 1) Each row in each data frame is a child.
# 2) Our linking variable stfcid is similarly named and encoded.
# 3) No other variable names match. 

# Now lets link by joining! 
# Although R will try to guess your linking variable, 
# it's always prudent to specify it explicitly. 
dlink <- ac_link |> 
  left_join(nc_link, 
            join_by(stfcid))

# Note a few things: 
nrow(dlink) == nrow(ac_link) 
ncol(dlink) == ncol(ac_link) + ncol(nc_link) - 1 # one linking variable

# We have a 7% failed link rate, 
# or false negative rate (i.e. is.na(nrep))
dlink |> 
  count(nrep) |> 
  mutate(pct = n/sum(n)*100)

# Note that we might also be concerned about report timing. 
# Because we include all maltreatment reports across FY2023, 
# for some children, later reports may represent 
# subsequent instances of maltreatment, and should be 
# subtracted out from our maltreatment history measures.  
dlink |>  
  filter(maxdt > dodfcdt) |> # report follows foster care discharge
  select(stfcid, maxdt, dodfcdt) |> 
  nrow()

# Nevertheless, we can now explore our research question
# What are the maltreatment histories of children 
# aged 0-3 placed into foster care in FY2023?
dlink |> 
  pivot_longer(cols = c(nrep, nsub), names_to = 'type', values_to = 'n') |> 
  ggplot(aes(x = n)) + 
  geom_histogram(binwidth = 1) + 
  facet_wrap(~type) + 
  scale_x_continuous(breaks = 0:10)

# Finally, let's save our data
write_rds(dlink,paste0(data,'linked_data.rds') )