lepidoptera is a domain-agnostic causal simulation engine for structural causal models with mixed variable types. The package separates the causal graph (which variables exist and who depends on whom) from the mechanism (how each variable is generated given its parents), enabling clean intervention semantics and individual-level counterfactual queries.
The core workflow has four steps:
feature_spec()).mechanism_library_get() or a custom factory).new_system(),
add_feature(), add_dependency(),
finalize_system()).simulate()) and propagate interventions
(propagate()).library("lepidoptera")
#>
#> Attaching package: 'lepidoptera'
#> The following object is masked from 'package:stats':
#>
#> simulate
# 1. Feature specs
age_spec <- feature_spec("age", "continuous", list(min = 18, max = 100))
cvd_spec <- feature_spec("cvd", "categorical", list(levels = c("no", "yes")))
# 2. Mechanisms
mech_age <- mechanism_library_get("normal")(mean = 55, sd = 12)
mech_cvd <- mechanism_library_get("multinomial_logit")(
levels = c("no", "yes"),
coefs = list(yes = list(age = 0.04)),
intercepts = list(yes = -3.0)
)
# 3. System
sys <- new_system() |>
add_feature(age_spec) |>
add_feature(cvd_spec) |>
add_dependency(character(0), "age", mech_age) |>
add_dependency("age", "cvd", mech_cvd) |>
finalize_system()
# 4. Simulate
pop <- simulate(sys, n = 1000, seed = 1)
head(pop)
#> age cvd
#> 1 61.22352 yes
#> 2 48.69175 no
#> 3 45.66221 no
#> 4 47.61104 no
#> 5 43.97003 no
#> 6 50.61631 yes
prop.table(table(pop$cvd))
#>
#> no yes
#> 0.689 0.311propagate() applies one or more interventions and
returns both baseline and intervened populations. With
mode = "full" each unit is drawn independently from the
interventional distribution (suitable for ATE estimation); with
mode = "partial" unaffected nodes retain their baseline
draws so the same individual appears in both worlds (suitable for ITE
estimation).
# Population-level intervention: what if everyone were 70?
result_full <- propagate(
sys,
interventions = list(do_fix("age", 70)),
n = 1000, seed = 1, mode = "full"
)
cat("Baseline CVD rate: ", mean(result_full$baseline$cvd == "yes"), "\n")
#> Baseline CVD rate: 0.311
cat("Under age=70 CVD rate: ", mean(result_full$intervened$cvd == "yes"), "\n")
#> Under age=70 CVD rate: 0.439
# Individual-level counterfactual: same 1000 people, age set to 70
result_cf <- propagate(
sys,
interventions = list(do_fix("age", 70)),
n = 1000, seed = 1, mode = "partial"
)
ite <- (result_cf$intervened$cvd == "yes") - (result_cf$baseline$cvd == "yes")
cat("Mean ITE: ", mean(ite), "\n")
#> Mean ITE: 0.128When observed data are available, copula_mechanism()
fits the full joint dependence structure non-parametrically and uses it
as the conditional generating rule at simulation time.
set.seed(42)
obs <- data.frame(
sbp = rnorm(300, 130, 15),
dbp = NA_real_
)
obs$dbp <- 0.5 * obs$sbp + rnorm(300, 10, 8)
sbp_spec <- feature_spec("sbp", "continuous", list(min = 70, max = 250))
dbp_spec <- feature_spec("dbp", "continuous", list(min = 40, max = 150))
mech_sbp <- mechanism_library_get("normal")(mean = 130, sd = 15)
mech_dbp <- copula_mechanism(
variables = c("dbp", "sbp"),
family = "auto",
fit_data = obs
)
sys_bp <- new_system() |>
add_feature(sbp_spec) |> add_feature(dbp_spec) |>
add_dependency(character(0), "sbp", mech_sbp) |>
add_dependency("sbp", "dbp", mech_dbp) |>
finalize_system()
sim_bp <- simulate(sys_bp, n = 500, seed = 7)
cat(sprintf(
"Observed Spearman rho: %.2f | Simulated: %.2f\n",
cor(obs$sbp, obs$dbp, method = "spearman"),
cor(sim_bp$sbp, sim_bp$dbp, method = "spearman")
))
#> Observed Spearman rho: 0.67 | Simulated: 0.62The companion paper (Sarti, 2026, Journal of Statistical Software) provides a full technical description of the framework, the mechanism library, the partial propagation algorithm, and three worked illustrations spanning clinical trial simulation, emergency-ward digital twins, and epidemiological mediation analysis.