Graph-derived piecewise imputation for missing data (drm

Missing data is pervasive in ecology, evolutionary biology, and the environmental sciences (Little and Rubin 2019; Enders 2010). In mediation chains and structural equation models ( $X \to M \to Y$ ), intermediate traits $M$ (such as hormone titers, metabolic rates, behavioural assays, or tissue samples) are often missing on a substantial fraction of individuals.

When mediator $M$ is missing conditionally on the outcome $Y$ (Missing at Random, MAR), dropping incomplete rows via complete-case analysis induces severe bias in both the direct and indirect paths.

drmSEM provides graph-derived piecewise imputation via impute = "auto" in drm_sem() and drm_psem(). This vignette explains:

Why complete-case analysis fails in causal mediation chains.
How drmSEM derives imputation models directly from the causal DAG without ad-hoc user choices.
Supported response and predictor families (Gamma, lognormal, student, nbinom2, ZIP, gaussian).
Multi-parent imputation ( $k = 2$ incomplete parents) and engine constraints.
Inspection tools: imputation(), imputed(), and check_sem().
The honest scientific boundary: why piecewise graph-derived imputation is not joint FIML, and how within-node uncertainty is propagated.

1. The problem: mediator missingness in causal chains

Consider a classic mediation chain where an environmental driver $X$ influences a physiological mediator $M$ , which in turn affects fitness $Y$ :

$X \to M \to Y$

# Simulate a mediation chain with outcome-dependent missingness (MAR)
n <- 400
x <- rnorm(n, mean = 0, sd = 1)
m_true <- 0.7 * x + rnorm(n, sd = 0.5)
y <- 0.5 * m_true + 0.4 * x + rnorm(n, sd = 0.5)

dat <- data.frame(x = x, m = m_true, y = y)

# Induce MAR missingness: individuals with low fitness (Y) are less likely to be sampled for M
p_miss <- plogis(-0.5 - 1.2 * scale(y))
dat$m[runif(n) < p_miss] <- NA

# Check missingness proportion
mean(is.na(dat$m))
#> [1] 0.38

If we fit this system with standard complete-case analysis (impute = "none"), any row with missing $M$ is dropped from the regression of $Y$ on $M$ and $X$ . Because missingness depends on $Y$ , complete-case estimates of $M \to Y$ and $X \to Y$ are substantially biased.

2. Graph-derived automatic imputation (`impute = "auto"`)

In traditional structural equation modelling, handling missing predictors requires either: - Specifying an ad-hoc imputation model (risking uncongeniality with the analysis model), or - Running external chained equations (MICE) with manual variable selections.

In drmSEM, the causal DAG already specifies the imputation model. Because $M$ is an endogenous node in the SEM with its own formula and family ( $M \sim X$ , $\text{family} = \dots$ ), drmSEM automatically borrows node $M$ ’s model and wires it directly into node $Y$ ’s likelihood as a modelled missing predictor (mi(m)):

sem_imputed <- drm_sem(
  m = drm_node(drmTMB::bf(m ~ x), family = stats::gaussian()),
  y = drm_node(drmTMB::bf(y ~ m + x), family = stats::gaussian()),
  data = dat,
  impute = "auto"
)
#> ℹ Imputation models derived from the graph: "y <- m".
#> ℹ Each is the incomplete parent's own node model. The downstream node
#>   re-estimates it rather than sharing the parent's estimates, so this is not
#>   full-information maximum likelihood across the SEM. See
#>   docs/design/13-missing-data.md.
#> Warning: Nodes would be fitted on different row sets.
#> ! Rows used of 400: m = 248, y = 400.
#> ℹ Imputed nodes keep their rows; a node whose own RESPONSE is incomplete still
#>   loses them, which is usually the difference you see above.
#> ℹ See `attr(x, "alignment_issues")`. Use `na_action = "common"` to fit every
#>   node on the shared complete-case set, or `"fail"` to make this an error.
#> ℹ Fitting node "m"
#> ✔ Fitting node "y" [104ms]
#> 
#> ℹ Fitting node "y"
#> ✔ Fitting node "y" [178ms]

summary(sem_imputed)
#> 
#> ── drmSEM summary ──────────────────────────────────────────────────────────────
#> 2 endogenous nodes; order "m" and "y"
#> 
#> ── Edges by component ──
#> 
#> mu: 3
#> 
#> ── Paths ──
#> 
#> <drmSEM paths: 3 component-labelled coefficients>
#>  from to component     link  term estimate std.error statistic  p.value
#>     x  m        mu identity     x   0.6308    0.0338    18.672 8.35e-78
#>     m  y        mu identity mi(m)   0.5026    0.0624     8.053 8.07e-16
#>     x  y        mu identity     x   0.4218    0.0508     8.300 1.05e-16
#>  endogenous
#>       FALSE
#>        TRUE
#>       FALSE

No manual imputation formulas or custom missing-data scripts are required: the DAG supplies the exact parent set and distributional family for $M$ .

3. Supported continuous and discrete families

drmSEM supports graph-derived imputation across a wide range of response and predictor distributions:

Continuous families

Gaussian (gaussian): Standard continuous location-scale responses.
Gamma (Gamma / gamma): Positive continuous traits (e.g. body size, metabolic rate, dispersal distance).
Log-normal (lognormal): Skewed positive measurements.
Student- $t$ (student): Heavy-tailed robust continuous outcomes.
Beta (beta): Proportions and bounded indices on $(0, 1)$ .

Discrete and count families

Poisson (poisson): Unbounded count data (e.g. egg counts, territory visits).
Negative Binomial (nbinom2): Overdispersed ecological count responses with quadratic variance $\mu + \mu^2/\theta$ .
Zero-Inflated Poisson (ZIP / zi_poisson): Count data with excess structural zeros (e.g. parasite counts, mating successes).
Binomial & Beta-Binomial (binomial, beta_binomial): Binary trials and overdispersed success proportions.

# Example with a Gaussian mediator and Negative Binomial count response
dat_count <- data.frame(
  temp = rnorm(n),
  x_env = rnorm(n)
)
dat_count$biomass <- 0.6 * dat_count$temp + rnorm(n, sd = 0.5)
dat_count$parasites <- rpois(n, lambda = exp(0.5 + 0.4 * dat_count$biomass - 0.3 * dat_count$x_env))

# Introduce missingness in biomass
dat_count$biomass[runif(n) < 0.25] <- NA

sem_count <- drm_sem(
  biomass   = drm_node(drmTMB::bf(biomass ~ temp), family = stats::gaussian()),
  parasites = drm_node(drmTMB::bf(parasites ~ biomass + x_env), family = drmTMB::nbinom2()),
  data = dat_count,
  impute = "auto"
)
#> ℹ Imputation models derived from the graph: "parasites <- biomass".
#> ℹ Each is the incomplete parent's own node model. The downstream node
#>   re-estimates it rather than sharing the parent's estimates, so this is not
#>   full-information maximum likelihood across the SEM. See
#>   docs/design/13-missing-data.md.
#> Warning: Nodes would be fitted on different row sets.
#> ! Rows used of 400: biomass = 300, parasites = 400.
#> ℹ Imputed nodes keep their rows; a node whose own RESPONSE is incomplete still
#>   loses them, which is usually the difference you see above.
#> ℹ See `attr(x, "alignment_issues")`. Use `na_action = "common"` to fit every
#>   node on the shared complete-case set, or `"fail"` to make this an error.
#> ℹ Fitting node "biomass"
#> ✔ Fitting node "parasites" [32ms]
#> 
#> ℹ Fitting node "parasites"
#> ✔ Fitting node "parasites" [568ms]

paths(sem_count)
#> <drmSEM paths: 3 component-labelled coefficients>
#>     from        to component     link        term estimate std.error statistic
#>     temp   biomass        mu identity        temp   0.6029    0.0281    21.422
#>  biomass parasites        mu      log mi(biomass)   0.4625    0.0514     9.000
#>    x_env parasites        mu      log       x_env  -0.2947    0.0405    -7.276
#>    p.value endogenous
#>  8.29e-102      FALSE
#>   2.26e-19       TRUE
#>   3.43e-13      FALSE

4. Multi-parent imputation ( $k = 2$ ) and engine constraints

When a downstream node $Y$ has two incomplete endogenous parents ( $X \to M_1 \to Y$ and $X \to M_2 \to Y$ ), drmSEM automatically derives independent modelled imputation terms for both parents:

# Two incomplete Gaussian parents
dat2 <- data.frame(x = rnorm(n))
dat2$m1 <- 0.6 * dat2$x + rnorm(n, sd = 0.5)
dat2$m2 <- -0.4 * dat2$x + rnorm(n, sd = 0.5)
dat2$y  <- 0.5 * dat2$m1 + 0.3 * dat2$m2 + 0.2 * dat2$x + rnorm(n, sd = 0.5)

# Missingness in both parents
dat2$m1[runif(n) < 0.20] <- NA
dat2$m2[runif(n) < 0.20] <- NA

sem_2p <- drm_sem(
  m1 = drm_node(drmTMB::bf(m1 ~ x), family = stats::gaussian()),
  m2 = drm_node(drmTMB::bf(m2 ~ x), family = stats::gaussian()),
  y  = drm_node(drmTMB::bf(y ~ m1 + m2 + x), family = stats::gaussian()),
  data = dat2,
  impute = "auto"
)
#> ℹ Imputation models derived from the graph: "y <- m1, m2".
#> ℹ Each is the incomplete parent's own node model. The downstream node
#>   re-estimates it rather than sharing the parent's estimates, so this is not
#>   full-information maximum likelihood across the SEM. See
#>   docs/design/13-missing-data.md.
#> Warning: Nodes would be fitted on different row sets.
#> ! Rows used of 400: m1 = 321, m2 = 320, y = 400.
#> ℹ Imputed nodes keep their rows; a node whose own RESPONSE is incomplete still
#>   loses them, which is usually the difference you see above.
#> ℹ See `attr(x, "alignment_issues")`. Use `na_action = "common"` to fit every
#>   node on the shared complete-case set, or `"fail"` to make this an error.
#> ℹ Fitting node "m1"
#> ✔ Fitting node "m2" [32ms]
#> 
#> ℹ Fitting node "m2"
#> ✔ Fitting node "y" [41ms]
#> 
#> ℹ Fitting node "y"
#> ✔ Fitting node "y" [157ms]

paths(sem_2p)
#> <drmSEM paths: 5 component-labelled coefficients>
#>  from to component     link   term estimate std.error statistic  p.value
#>     x m1        mu identity      x   0.5426    0.0277    19.586 2.04e-85
#>     x m2        mu identity      x  -0.3609    0.0293   -12.308 8.21e-35
#>    m1  y        mu identity mi(m1)   0.4636    0.0571     8.121 4.61e-16
#>    m2  y        mu identity mi(m2)   0.1805    0.0564     3.202 1.36e-03
#>     x  y        mu identity      x   0.2247    0.0464     4.847 1.26e-06
#>  endogenous
#>       FALSE
#>       FALSE
#>        TRUE
#>        TRUE
#>       FALSE

Engine constraints and clear error reporting

drmSEM enforces strict honesty regarding engine limits: - $k > 2$ incomplete parents: The fitting engine currently limits simultaneous modelled imputation to $k \le 2$ terms per node. Attempting $k > 2$ aborts with an actionable error directing the user to complete or impute upstream predictors. - $k = 2$ on non-Gaussian responses: Multi-parent imputation currently requires a Gaussian response node; non-Gaussian $k = 2$ cells fail loudly with clear diagnostics.

5. Inspecting imputation models and diagnostics

drmSEM provides three dedicated diagnostic accessors:

1. `imputation(sem)`

Summarizes the exact graph-derived imputation plan for each node:

imputation(sem_imputed)
#>   node variable model   family n_missing uncertainty_status std_error_usable
#> 1    y        m m ~ x gaussian       152                 ok             TRUE

2. `imputed(sem)`

Extracts observation-level imputed predictions, standard errors, and status flags:

imp_vals <- imputed(sem_imputed)
head(imp_vals)
#>   node variable original_row model_row observed   estimate std_error
#> 1    y        m            2         2    FALSE  1.2540664 0.4458052
#> 2    y        m            4         4    FALSE -1.2781684 0.4531732
#> 3    y        m            5         5    FALSE -0.5646132 0.4470518
#> 4    y        m            7         7    FALSE -0.8026063 0.4493321
#> 5    y        m            8         8    FALSE -0.8395815 0.4467388
#> 6    y        m            9         9    FALSE -0.6158088 0.4465449
#>             source uncertainty_status
#> 1 conditional_mode                 ok
#> 2 conditional_mode                 ok
#> 3 conditional_mode                 ok
#> 4 conditional_mode                 ok
#> 5 conditional_mode                 ok
#> 6 conditional_mode                 ok

Observations that were originally missing carry observed = FALSE and report imputed values along with their model-based standard errors.

3. `check_sem(sem)`

Performs a comprehensive health check across the entire SEM, reporting sample sizes, convergence status, covariance availability, and missingness:

check_sem(sem_imputed)
#> 
#> ── drmSEM diagnostics ──
#> 
#>  node   family components nobs converged vcov_available sampler
#>     m gaussian         mu  248      TRUE           TRUE    TRUE
#>     y gaussian         mu  400      TRUE           TRUE    TRUE
#> Exogenous variables: "x"
#> Warning: Nodes were fitted on different numbers of observations.
#> ℹ Path coefficients then come from different samples. Refit with `na_action =
#>   "common"` in `drm_sem()` to use one shared complete-case set.

6. Honest piecewise bounds versus FIML

It is critical to distinguish piecewise graph-derived imputation from Full-Information Maximum Likelihood (FIML):

Within-node parameter estimation: In piecewise SEM, node $Y$ re-estimates the missing parent model ( $M \sim X$ ) within its own likelihood function. It does not rigidly fix the parameters to node $M$ ’s point estimates.
Uncertainty propagation: Imputation uncertainty for $M$ is fully accounted for within node $Y$ ’s joint Hessian (standard errors for $M \to Y$ and $X \to Y$ correctly reflect missing-data uncertainty).
No across-node joint likelihood: Because drmSEM is piecewise, it does not fit a single global joint covariance matrix across all endogenous nodes. It must never be described as whole-SEM FIML.

This approach provides a robust, principled solution for missing mediators that eliminates MAR attrition bias while respecting the piecewise architecture.

Summary of best practices

Activate auto-imputation: Use drm_sem(..., impute = "auto") whenever endogenous mediators contain missing observations.
Check the imputation plan: Run imputation(sem) to verify that derived parent formulas match your theoretical DAG.
Inspect individual imputations: Use imputed(sem) to verify that imputed values lie within reasonable biological bounds.
Audit model health: Run check_sem(sem) to confirm node convergence and Hessian positive-definiteness across all sub-models.

References

Enders, Craig K. 2010. Applied Missing Data Analysis. Guilford Press.

Little, Roderick J. A., and Donald B. Rubin. 2019. Statistical Analysis with Missing Data. 3rd ed. John Wiley & Sons.

Graph-derived piecewise imputation for missing data (drm_sem)