Skip to contents

Getting started with MORIE in R

Source: vignettes/intro.Rmd
intro.Rmd

Overview

rmorie is the R package of the MORIE toolkit. It covers dataset loading, causal estimators, survey sampling and weighting, signal processing, and the MRM (Multilevel Reconciliation Methodology) helpers for Canadian carceral, police and oversight data.

This vignette walks through a minimal session: load the package, look at the dataset catalogue, load one dataset, and run an average-treatment-effect estimator on a small synthetic example. The mrm-otis-walkthrough vignette runs the ten-estimator OTIS ensemble of morie_otis_analyze_a01_mrm().

Loading the package

The dataset catalogue

morie_dataset_catalog() returns a data frame with one row per dataset that morie_load_dataset() can load by key: Statistics Canada and Health Canada PUMFs, OTIS tables, Toronto Police layers and others. Reading the catalogue needs no network access.

catalog <- morie_dataset_catalog()
nrow(catalog)
#> [1] 71
table(catalog$source)
#> 
#>       cihi        hib       naps         oc       otis rmoriedata    statcan 
#>          6         14         24         12          5          2          1 
#>        tps        vsr 
#>          3          4
head(catalog[, c("key", "name", "source", "year")])
#>       key                      name source      year
#> 1   ocp21      CPADS 2021-2022 PUMF     oc 2021-2022
#> 2   occ22        CCS 2018-2022 PUMF     oc 2018-2022
#> 3   occ23             CCS 2023 PUMF     oc      2023
#> 4   occ24             CCS 2024 PUMF     oc      2024
#> 5 ocs22mf      CSADS 2021-2022 PUMF     oc 2021-2022
#> 6 ocs22bt CSADS 2021-2022 Bootstrap     oc 2021-2022

morie_dataset_info() returns the catalogue entry for one key as a list, including where the data are fetched from:

info <- morie_dataset_info("ocp21")
info[c("key", "name", "survey", "year", "ckan_resource_id")]
#> $key
#> [1] "ocp21"
#> 
#> $name
#> [1] "CPADS 2021-2022 PUMF"
#> 
#> $survey
#> [1] "cpads"
#> 
#> $year
#> [1] "2021-2022"
#> 
#> $ckan_resource_id
#> [1] "d2639429-c304-45a6-90b3-770562f4d46d"

Loading a dataset

morie_load_dataset() returns a data frame for any key in the catalogue. Remote datasets are downloaded on the first call and kept in a per-user cache under morie_cache_dir(), so later calls read the local copy. The chunk below is not run when the vignette is built, because the first call downloads the CPADS 2021–2022 PUMF (about 40 MB) from open.canada.ca:

df <- morie_load_dataset("ocp21")
dim(df)
# 40931   394

A simple ATE estimate

For users who already have a treatment / outcome / covariate dataset in hand, the estimators are designed to work on any tibble or data frame — there is no hard-coded column-name convention. The example below is fully synthetic and runnable without any external data.

set.seed(2026)
n <- 500
X1 <- rnorm(n)
X2 <- rnorm(n)
# Confounded treatment assignment.
treat <- as.integer(plogis(0.5 * X1 - 0.3 * X2) > runif(n))
# Outcome with a true ATE of +1.0 plus covariate effects.
y <- 1.0 * treat + 0.7 * X1 - 0.2 * X2 + rnorm(n, sd = 0.5)

df_synth <- data.frame(y = y, treat = treat, X1 = X1, X2 = X2)
result <- morie_estimate_ate(
  data       = df_synth,
  outcome    = "y",
  treatment  = "treat",
  covariates = c("X1", "X2")
)
print(result)
#> $n_used
#> [1] 500
#> 
#> $n_discarded
#> [1] 0
#> 
#> $estimand
#> [1] "ATE on the full sample"
#> 
#> $ate
#> [1] 0.9465714
#> 
#> $se
#> [1] 0.08977947
#> 
#> $ci_lower
#> [1] 0.7706036
#> 
#> $ci_upper
#> [1] 1.122539
#> 
#> $n
#> [1] 500
#> 
#> $ess
#> [1] 461.3557

The returned object is a list: the point estimate ate, its standard error se, the confidence limits ci_lower and ci_upper, the number of rows used and discarded by propensity trimming, and the effective sample size of the weights (ess).

Companion estimators

morie_estimate_att(), morie_estimate_atc(), and morie_estimate_aipw() follow the same calling convention. The augmented IPW estimator (morie_estimate_aipw()) is consistent if either the propensity model or the outcome model is correctly specified.

result_aipw <- morie_estimate_aipw(
  data       = df_synth,
  outcome    = "y",
  treatment  = "treat",
  covariates = c("X1", "X2")
)
print(result_aipw)
#> $n_used
#> [1] 500
#> 
#> $n_discarded
#> [1] 0
#> 
#> $estimand
#> [1] "ATE on the full sample"
#> 
#> $ate
#> [1] 0.9494884
#> 
#> $se
#> [1] 0.0451242
#> 
#> $ci_lower
#> [1] 0.8610449
#> 
#> $ci_upper
#> [1] 1.037932
#> 
#> $n
#> [1] 500

Where to go next