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, 
+                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) 
> # 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'))
> # 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 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()
             [,1]
subyr        0.00
staterr      0.00
stfcid       0.00
afcarsid     0.00
rptid        0.00
chid         0.00
rptdt        0.00
rptdisp      0.11
rpdispdt     0.00
chprior      0.00
fcmoney     88.69
per1rel     73.55
per1rel_new 73.85
> # 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
    per1rel      n
      <int>  <int>
 1:       1 158461
 2:       2   5924
 3:       3    248
 4:       4    370
 5:       5    923
 6:       6    830
 7:       7  11938
 8:       8   1407
 9:       9    714
10:      10    484
11:      33     18
12:      88   6428
13:      99   2202
14:      NA 528091
> nc |> 
+   count(per1rel_new) # clean encoding
                   per1rel_new      n
                        <char>  <int>
1:               Foster parent    636
2:             Friend/neighbor    484
3:              Legal guardian   1407
4:                       Other   6428
5:  Other relative, non-foster   5924
6:                      Parent 158461
7:                Professional   2467
8: Unmarried partner of parent  11938
9:                        <NA> 530293
> # 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))
Warning message:
In geom_histogram(stat = "identity") :
  Ignoring unknown parameters: `binwidth` and `bins`

> # 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()
            [,1]
fy          0.00
st          0.00
recnumbr    0.00
dob         0.00
totalrem    0.00
rem1dt      0.00
latremdt    0.00
ageatlatrem 0.00
dodfcdt     0.32
entered     0.00
exited      0.00
inatend     0.00
inatstart   0.00
stfcid      0.00
version     0.37
nsub        0.03
nrep        0.03
maxdt       0.03
> # 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)
+          )
> # 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)

Call:
glm(formula = subind ~ chprior + fcmoney, family = "binomial", 
    data = ncs)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  0.10250    0.03513   2.918  0.00353 ** 
chpriorNo    0.03731    0.01857   2.009  0.04451 *  
fcmoneyNo   -0.59147    0.03475 -17.020  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 75424  on 56193  degrees of freedom
Residual deviance: 75133  on 56191  degrees of freedom
  (522063 observations deleted due to missingness)
AIC: 75139

Number of Fisher Scoring iterations: 4

> nobs(m_cc)
[1] 56194
> nrow(ncs)
[1] 578257
> # 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) 
        subyr staterr         stfcid     afcarsid        rptid
        <int>  <char>         <char>        <i64>       <char>
     1:  2022      CT CT000410019186    410019186 100005B45C09
     2:  2021      CT CT000410056939    410056939 1000057AF159
     3:  2020      CT CT000410059634    410059634 1000053A36B1
     4:  2020      CT CT000410072065    410072065 1000053A3338
     5:  2021      CT CT000410513558    410513558 1000057AC1F6
    ---                                                       
578253:  2020      VT VT999822899999 999822899999 1000053D719F
578254:  2021      VT VT999822899999 999822899999 10000568BCDB
578255:  2020      VT VT999912899999 999912899999 1000053D7766
578256:  2022      VT VT999912899999 999912899999 100005AF3D3D
578257:  2020      VT VT999922899999 999922899999 1000053D81A7
                chid      rptdt         rptdisp   rpdispdt chprior
              <char>     <IDat>          <fctr>     <IDat>  <fctr>
     1: 100001862356 2022-04-08   Substantiated 2022-04-08     Yes
     2: 100001816BC8 2021-01-23   Substantiated 2021-02-14     Yes
     3: 10000184E831 2020-01-23 Unsubstantiated 2020-03-04     Yes
     4: 10000539F44A 2020-02-08 Unsubstantiated 2020-03-06      No
     5: 1000057A811C 2020-12-23 Unsubstantiated 2021-01-15      No
    ---                                                           
578253: 1000053D6607 2019-12-23 Unsubstantiated 2020-01-27      No
578254: 1000053D6607 2021-09-08 Unsubstantiated 2021-09-17      No
578255: 1000053D6609 2020-03-08   Substantiated 2020-03-31      No
578256: 1000053D6609 2022-05-23 Unsubstantiated 2022-06-02     Yes
578257: 1000053D660A 2019-12-08 Unsubstantiated 2020-03-04      No
        fcmoney per1rel per1rel_new subind
         <fctr>   <int>      <char>  <num>
     1:    <NA>       1      Parent      1
     2:    <NA>      NA        <NA>      1
     3:    <NA>      NA        <NA>      0
     4:    <NA>      NA        <NA>      0
     5:    <NA>      NA        <NA>      0
    ---                                   
578253:      No      NA        <NA>      0
578254:      No      NA        <NA>      0
578255:      No       1      Parent      1
578256:      No      NA        <NA>      0
578257:      No      NA        <NA>      0
> # 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')
   subyr staterr         stfcid     afcarsid        rptid         chid
   <int>  <char>         <char>        <i64>       <char>       <char>
1:  2022      VT VT582102899999 582102899999 100005AF3B72 10000568A272
        rptdt         rptdisp   rpdispdt chprior fcmoney per1rel
       <IDat>          <fctr>     <IDat>  <fctr>  <fctr>   <int>
1: 2022-01-23 Unsubstantiated 2022-03-01     Yes    <NA>      NA
   per1rel_new subind
        <char>  <num>
1:        <NA>      0
> # 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()
> # 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)
Warning message:
Number of logged events: 1 

> # 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)
Class: mids
Number of multiple imputations:  5 
Imputation methods:
   subyr  staterr    rptdt   stfcid   subind  chprior  fcmoney 
      ""       ""       ""       ""    "pmm"       "" "logreg" 
PredictorMatrix:
        subyr staterr rptdt stfcid subind chprior fcmoney
subyr       0       0     0      0      0       0       0
staterr     0       0     0      0      0       0       0
rptdt       0       0     0      0      0       0       0
stfcid      0       0     0      0      0       0       0
subind      0       0     0      0      0       1       0
chprior     0       0     0      0      0       0       0
Number of logged events:  1 
  it im dep     meth     out
1  0  0     constant staterr
> # We can easily diagnose the convergence of our imputation model 
> # (see additional resources).  
> plot(imp)
> # 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 %`)
Waiting for profiling to be done...
> 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')
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.