Skip to contents

The native engines: DAG to publication table without leaving rmorie

Source: vignettes/native-engines.Rmd
native-engines.Rmd

This vignette walks a whole causal pipeline on one simulated study using only library(rmorie). No MatchIt, no survey, no DoubleML, no dagitty, no fixest: each step below runs on an engine in this package. The source repository’s tests/cross/ directory holds the tests that compare these engines with the reference packages (for example test-morie_vs_matchit.R, test-morie_vs_dagitty.R, test-morie_vs_doubleml.R).

The data

A simulated observational study: three confounders drive both a binary treatment and a continuous outcome, with a true effect of 0.9.

n <- 600
age  <- rnorm(n); educ <- rnorm(n); re74 <- rnorm(n)
treat <- rbinom(n, 1, plogis(0.4 * age + 0.3 * educ + 0.3 * re74))
re78  <- 1 + 0.9 * treat + 0.5 * age + 0.4 * educ + 0.6 * re74 +
  rnorm(n)
d <- data.frame(age, educ, re74, treat, re78)

0. Guard the categories first (module 25)

Before any model: audit the categorical columns. Silent category-mapping errors — numeric codes imported without their labels, positional recodes, alphabetical releveling — can corrupt a disparity analysis without any error message. rmorie’s audit flags such columns, morie_safe_recode() recodes by label rather than by position, and morie_crosstab_verify() checks a recode against its mapping (it stops on a mismatch):

grp_codes <- sample(c("1", "2"), n, replace = TRUE)   # an import hazard
morie_audit_categories(data.frame(group = factor(grp_codes)))
#> Categorical audit: 1 column(s)
#>   group            factor     2 level(s), reference '1'
#>     !! HAZARD: all labels numeric-looking (1,2...): likely imported CODES whose value labels were lost; as.numeric() on this column returns level INDICES, not data

grp <- morie_safe_recode(grp_codes, c("1" = "Control", "2" = "Program"))
morie_crosstab_verify(grp_codes, grp,
                      c("1" = "Control", "2" = "Program"))

1. The DAG (module 13) 

dag <- morie_dag(
  edges = c("re74 -> treat", "re74 -> re78",
            "age -> treat",  "age -> re78",
            "educ -> treat", "educ -> re78",
            "treat -> re78"),
  exposure = "treat", outcome = "re78")
dag
#> morie causal DAG: 5 nodes, 7 edges
#>   exposure: treat  outcome: re78 
#>    re74 -> treat 
#>    re74 -> re78 
#>    age -> treat 
#>    age -> re78 
#>    educ -> treat 
#>    educ -> re78 
#>    treat -> re78

2. Identification

Backdoor identification runs Bayes-Ball d-separation with the canonical adjustment set (validated against dagitty in tests/cross/):

id <- morie_dag_identify(dag)
id$identified
#> [1] TRUE
id$adjustment_set
#> [1] "re74" "age"  "educ"

3. Matching on the identified set (module 1) 

m <- morie_matching_nearest_neighbor(d, "treat", id$adjustment_set,
                                     replace = TRUE)
#> Warning: Fewer control units than treated units; not all treated units will get
#> a match.
nrow(m$match_pairs)
#> [1] 304
morie_matching_balance(m$matched_data, "treat",
                       id$adjustment_set)$max_smd
#> [1] 0.2750579

4. Estimation through the DAG (modules 10 + 13)

morie_dag_estimate() routes the identified set into the native estimators — here double machine learning (cross-fit ridge + logit nuisances, Neyman-orthogonal score):

ate <- morie_dag_estimate(dag, d, method = "backdoor.dml")
unlist(ate[c("ate", "se", "ci_lower", "ci_upper")])
#>        ate         se   ci_lower   ci_upper 
#> 0.91198904 0.08765135 0.74019555 1.08378253

5. Refutation (module 13)

DoWhy-style robustness checks, all native — a placebo treatment should erase the effect; subsets and a random common cause should leave it stable:

for (mth in c("placebo_treatment", "random_common_cause",
              "data_subset")) {
  rf <- morie_dag_refute(dag, d, mth, n_reps = 10L)
  cat(sprintf("%-22s passed: %s\n", mth, rf$passed))
}
#> placebo_treatment      passed: TRUE
#> random_common_cause    passed: TRUE
#> data_subset            passed: TRUE

6. The MRM pipeline and a publication table (module 24)

The MRM flagship composes several native estimators on the same specification, applies a multiple-testing correction across them, pools an inverse-variance consensus, and renders the table with rmorie’s own renderer (text, markdown, LaTeX, or HTML):

eff <- morie_mrm_estimate_causal_effect(
  d, "treat", "re78", id$adjustment_set,
  methods = c("matching", "ate", "dml"))
#> Warning: Fewer control units than treated units; not all treated units will get
#> a match.
cat(morie_mrm_report(eff, format = "markdown"), sep = "\n")
#> | Method                   | Estimate | SE     | 95% CI        | p (adj.) |
#> |--------------------------|----------|--------|---------------|----------|
#> | matching (rmorie native) | 0.822*** | 0.142  | [0.543,  1.1] | 7.76e-09 |
#> | ipw ate (rmorie native)  | 0.865*** | 0.135  | [0.601, 1.13] | 2.91e-10 |
#> | dml plr (rmorie native)  | 0.912*** | 0.0877 | [0.74, 1.08]  | 7.08e-25 |
#> Consensus (inverse-variance): 0.882 (SE 0.11). holm correction; n = 600.

The consensus standard error is the inverse-variance-weighted mean of the estimators’ standard errors, not sqrt(1 / sum(w)): the estimators share one data set, so their errors are correlated and the independent-studies formula would be too small.

Where to look next