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. > ## 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) > # 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')) > # Recall from S2 that nc is already "stacked" > # after we unlisted our NCANDS submission-year > # into a single data frame. > nc |> + count(staterr,subyr) staterr subyr n 1: CT 2020 74880 2: CT 2021 79922 3: CT 2022 79203 4: CT 2023 76600 5: CT 2024 70930 6: MA 2020 16093 7: MA 2021 15555 8: MA 2022 17957 9: MA 2023 19654 10: MA 2024 20495 11: ME 2020 16165 12: ME 2021 14249 13: ME 2022 15501 14: ME 2023 15644 15: ME 2024 14980 16: NH 2020 25308 17: NH 2021 23445 18: NH 2022 20607 19: NH 2023 22437 20: NH 2024 20134 21: RI 2020 3622 22: RI 2021 3371 23: RI 2022 4574 24: RI 2023 4744 25: RI 2024 4207 26: VT 2020 9308 27: VT 2021 8078 28: VT 2022 7009 29: VT 2023 7215 30: VT 2024 6151 staterr subyr n > # 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) st fy n 1: CT 2023 12359 2: CT 2024 12147 3: MA 2023 4593 4: MA 2024 4687 5: ME 2023 1852 6: ME 2024 1786 7: NH 2023 3396 8: NH 2024 3442 9: RI 2023 1530 10: RI 2024 1428 11: VT 2023 2636 12: VT 2024 2412 > # 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)) + + geom_point() + + geom_line() > # 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 fy st stfcid entered exited inatend inatstart 1: 2023 CT CT000490093038 1 0 1 0 2: 2024 CT CT000490093038 0 0 1 1 3: 2023 CT CT000490273805 1 0 1 0 4: 2024 CT CT000490273805 0 1 0 1 5: 2023 CT CT000490376258 1 1 0 0 6: 2023 CT CT000490435166 1 0 1 0 7: 2024 CT CT000490435166 0 1 0 1 8: 2023 CT CT000490437443 0 0 1 1 9: 2024 CT CT000490437443 0 1 0 1 10: 2023 CT CT000490850584 0 0 1 1 > # 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))) + # visualize + geom_point() + + scale_x_continuous(limits = c(0,1)) + } > ac |> + link_success() `summarise()` has regrouped the output. ℹ Summaries were computed grouped by st, fy, and linked. ℹ Output is grouped by st and fy. ℹ Use `summarise(.groups = "drop_last")` to silence this message. ℹ Use `summarise(.by = c(st, fy, linked))` for per-operation grouping instead. > # 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() `summarise()` has regrouped the output. ℹ Summaries were computed grouped by st, fy, and linked. ℹ Output is grouped by st and fy. ℹ Use `summarise(.groups = "drop_last")` to silence this message. ℹ Use `summarise(.by = c(st, fy, linked))` for per-operation grouping instead. > # 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))) + # visualize + geom_point() + + scale_x_continuous(limits = c(.9,1)) `summarise()` has regrouped the output. ℹ Summaries were computed grouped by st, fy, and dob_match. ℹ Output is grouped by st and fy. ℹ Use `summarise(.groups = "drop_last")` to silence this message. ℹ Use `summarise(.by = c(st, fy, dob_match))` for per-operation grouping instead. > 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 Time difference of 7.028203 secs > #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) # A tibble: 18 × 2 nrep n 1 1 237029 2 2 66065 3 3 26751 4 4 12086 5 5 5795 6 6 2892 7 7 1276 8 8 666 9 9 338 10 10 165 11 11 74 12 12 36 13 13 13 14 14 12 15 15 3 16 16 2 17 17 1 18 13545 1 > nc_link |> count(nsub) # A tibble: 13 × 2 nsub n 1 0 209413 2 1 112666 3 2 22118 4 3 6241 5 4 1911 6 5 609 7 6 167 8 7 44 9 8 22 10 9 10 11 10 2 12 11 1 13 1822 1 > nc_link |> filter(nrep > 20 | nsub > 20) # A tibble: 1 × 4 stfcid nsub nrep maxdt 1 RI 1822 13545 2023-09-23 > head(nc_link) # A tibble: 6 × 4 stfcid nsub nrep maxdt 1 CT000410019186 1 1 2022-04-08 2 CT000410056939 1 1 2021-01-23 3 CT000410059634 0 1 2020-01-23 4 CT000410072065 0 1 2020-02-08 5 CT000410513558 0 1 2020-12-23 6 CT000410515973 1 1 2021-02-08 > # 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) fy st recnumbr dob totalrem rem1dt latremdt ageatlatrem 1: 2023 MA 1109152070 2020-12-15 2 2021-02-23 2023-03-07 2 2: 2023 MA 1109160599 2020-09-15 1 2022-11-28 2022-11-28 2 3: 2023 MA 1109208253 2020-09-15 1 2022-11-10 2022-11-10 2 4: 2023 MA 1109211241 2020-12-15 1 2023-05-05 2023-05-05 2 5: 2023 MA 1109230078 2021-01-15 1 2023-02-01 2023-02-01 2 6: 2023 MA 1109230986 2020-08-15 1 2023-04-10 2023-04-10 2 dodfcdt entered exited inatend inatstart stfcid version 1: 2023-08-22 1 1 0 0 MA001109152070 NA 2: 1 0 1 0 MA001109160599 NA 3: 1 0 1 0 MA001109208253 NA 4: 1 0 1 0 MA001109211241 NA 5: 1 0 1 0 MA001109230078 NA 6: 1 0 1 0 MA001109230986 NA > # 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) [1] TRUE > ncol(dlink) == ncol(ac_link) + ncol(nc_link) - 1 # one linking variable [1] TRUE > # 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) nrep n pct 1: 1 1266 47.08069914 2: 2 562 20.89996281 3: 3 334 12.42097434 4: 4 165 6.13611008 5: 5 86 3.19821495 6: 6 47 1.74786166 7: 7 23 0.85533656 8: 8 11 0.40907401 9: 9 2 0.07437709 10: 10 2 0.07437709 11: 11 2 0.07437709 12: NA 189 7.02863518 > # 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() [1] 51 > # 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) Warning message: Removed 378 rows containing non-finite outside the scale range (`stat_bin()`).