Skip to contents

MRM walkthrough on OTIS provincial data

Source: vignettes/mrm-otis-walkthrough.Rmd
mrm-otis-walkthrough.Rmd

What this vignette covers

The Offender Tracking Information System (OTIS) tables, published by the Ontario Ministry of the Solicitor General on data.ontario.ca, describe restrictive confinement and segregation in Ontario’s provincial institutions. This vignette shows the rmorie functions that work with them:

The OTIS analysis functions return lists with a tables element. This helper turns one of those tables into a data frame for printing:

as_table <- function(tab) {
  rows <- lapply(tab$rows, function(r) unlist(lapply(r, as.character)))
  out <- as.data.frame(do.call(rbind, rows), stringsAsFactors = FALSE)
  names(out) <- tab$headers
  out
}

Where the OTIS data come from

The dataset catalogue lists the OTIS tables that morie_load_dataset() fetches from data.ontario.ca:

cat_otis <- morie_dataset_catalog()
cat_otis[cat_otis$source == "otis", c("key", "name")]
#>        key
#> 37 otisloc
#> 38 otisa01
#> 39 otisb01
#> 40 otisb09
#> 41 otisc11
#>                                                                      name
#> 37 OTIS: Ontario's adult provincial correctional institutions (locations)
#> 38                   OTIS a01: Restrictive Confinement - Detailed Dataset
#> 39                               OTIS b01: Segregation - Detailed Dataset
#> 40            OTIS b09: Individuals in Segregation - Number of Placements
#> 41            OTIS c11: Individuals in Segregation/RC by Aggregate Length

The full tables are downloaded on the first call and then read from the local cache. This chunk is not run when the vignette is built: it needs the network, and the ensemble below takes a long time on the full a01 table (76,934 rows in October 2026; a test run had not finished within 20 minutes):

a01 <- morie_load_dataset("otisa01")
r_full <- morie_otis_analyze_a01_mrm(a01)

Offline, three kinds of OTIS data are available:

  • morie_sample("otis_b01") is a 1,000-row extract of table b01; morie_sample("otis_b09") and morie_sample("otis_c11") are the whole b09 and c11 tables.
  • morie_otis_load() returns the cached a01 table if there is one, and otherwise a 5-row synthetic sample in the a01 layout. Its numbers are not findings.
  • morie_synth_otis() simulates OTIS tables with the published column layout. The examples below use it for the a01 design.
b01 <- morie_sample("otis_b01")
dim(b01)
#> [1] 1000   18
a01_syn <- morie_synth_otis("a01", n = 300L, seed = 1L)
head(a01_syn[, c("EndFiscalYear", "UniqueIndividual_ID",
                 "MentalHealth_Alert", "SuicideRisk_Alert",
                 "SuicideWatch_Alert")], 3)
#>   EndFiscalYear UniqueIndividual_ID MentalHealth_Alert SuicideRisk_Alert
#> 1          2018           syn-00031                 No                No
#> 2          2021           syn-00012                 No                No
#> 3          2024           syn-00044                 No               Yes
#>   SuicideWatch_Alert
#> 1                 No
#> 2                Yes
#> 3                 No

The canonical a01 design

OTIS reassigns UniqueIndividual_ID every fiscal year, so the unit of analysis is the person-year. For each person-year, morie_otis_make_pair_alert_to_volatility_ruhela() builds:

  • ac: the number of distinct combinations of the three alert flags (mental health, suicide risk, suicide watch) across that person-year’s rows;
  • the treatment T_high_ac = 1 when ac >= 2;
  • the outcome Y_vm_count: the number of region changes, counting a change between the region at placement and the most recent region within a row, and a change of region from one row to the next;
  • the covariates Gender, Age_Category and EndFiscalYear.
pair <- morie_otis_make_pair_alert_to_volatility_ruhela(a01_syn)
nrow(pair$data)
#> [1] 227
table(pair$data$T_high_ac)
#> 
#>   0   1 
#> 174  53

The ten-estimator ensemble

morie_otis_analyze_a01_mrm() (an alias of morie_otis_analyze_a01_ruhela_formulations()) runs ten estimators on this design: IPW (Hajek), AIPW, g-computation, 1:1 propensity-score matching, propensity-score subclassification, IRM double machine learning (ATE, ATTE and ATC), matching followed by IRM-DML, an AIPW estimate of the ATC, the partially linear DML model, and a SuperLearner-stacked AIPW. Here it runs on the simulated a01 table, so the estimates describe the simulation, not Ontario:

res <- morie_otis_analyze_a01_mrm(a01_syn)
class(res)
#> [1] "morie_otis_analysis_result" "morie_rich_result"         
#> [3] "list"
str(res$summary_lines[c("n", "Treated prevalence")])
#> List of 2
#>  $ n                 : int 227
#>  $ Treated prevalence: chr "23.3%"
ens <- as_table(res$tables[[1]])
ens[, c("Estimator", "Estimand", "Estimate", "SE", "95% CI")]
#>                        Estimator Estimand Estimate     SE           95% CI
#> 1                 1. IPW (Hajek)      ATE  +2.0426 0.1572 [+1.735, +2.351]
#> 2               2. AIPW (RRZ DR)      ATE  +2.0494 0.1819 [+1.693, +2.406]
#> 3               3. g-computation      ATE  +2.0873 0.1472 [+1.799, +2.376]
#> 4                  4. PSM 1:1 NN      ATT  +2.1321 0.1525 [+1.833, +2.431]
#> 5            5. PSM subclass (5)      ATE  +2.0976 0.1651 [+1.774, +2.421]
#> 6                     6. IRM-DML      ATE  +2.0644 0.0903 [+1.887, +2.241]
#> 7                     6. IRM-DML     ATTE  +2.1322 0.1493 [+1.839, +2.425]
#> 8                     6. IRM-DML      ATC  +2.0438 0.1079 [+1.832, +2.255]
#> 9  7. PSM->IRM-DML (match_first)      ATE  +2.1679 0.1398 [+1.894, +2.442]
#> 10                 8. ATC (AIPW)      ATC  +2.0465 0.2128 [+1.629, +2.464]
#> 11    9. PLR (Chernozhukov 2018)    theta  +2.0852 0.1503 [+1.791, +2.380]
#> 12         10. SuperLearner AIPW      ATE  +2.0344 0.1528 [+1.735, +2.334]

The second table compares four standard errors for the IRM-DML estimate: pooled, clustered on fiscal year, clustered on person, and two-way:

as_table(res$tables[[2]])[, c("SE flavour", "ATE", "SE", "95% CI")]
#>                     SE flavour     ATE     SE           95% CI
#> 1                 pooled (iid) +2.0644 0.1801 [+1.711, +2.417]
#> 2       cluster: EndFiscalYear +2.0644 0.0903 [+1.887, +2.241]
#> 3 cluster: UniqueIndividual_ID +2.0644 0.1900 [+1.692, +2.437]
#> 4           multi-way: yr x id +2.0644 0.1093 [+1.850, +2.279]

The third table is a sensitivity arm built by morie_otis_make_pair_alert_to_volatility_naive(): the treatment is “at least two alert flags at once in some row of the person-year”, and the outcome is binary (any row whose two region columns differ):

as_table(res$tables[[3]])[, c("Estimator", "ATE", "SE", "95% CI")]
#>   Estimator     ATE     SE           95% CI
#> 1       IPW +0.0530 0.0504 [-0.046, +0.152]
#> 2      AIPW +0.0492 0.0539 [-0.056, +0.155]
#> 3    g-comp +0.0534 0.0533 [-0.051, +0.158]

The individual estimator outputs are kept in res$payload$ensemble.

The general multi-estimator entry point

morie_mrm_estimate_causal_effect() takes any data frame with a binary treatment and runs up to four estimators (nearest-neighbour matching, IPW, AIPW and DML) on the same specification. It adjusts the p-values for multiple testing (Holm by default) and pools the estimates. morie_mrm_report() renders the result:

eff <- morie_mrm_estimate_causal_effect(
  pair$data, treatment = pair$T, outcome = pair$Y,
  covariates = pair$covariates,
  methods = c("matching", "ate", "aipw", "dml"),
  correction = "holm"
)
cat(morie_mrm_report(eff, format = "markdown"), sep = "\n")
#> | Method                   | Estimate | SE    | 95% CI       | p (adj.) |
#> |--------------------------|----------|-------|--------------|----------|
#> | matching (rmorie native) | 2.06***  | 0.185 | [1.69, 2.42] | 1.22e-28 |
#> | ipw ate (rmorie native)  | 2.05***  | 0.159 | [1.74, 2.36] | 5.6e-38  |
#> | aipw (rmorie native)     | 2.06***  | 0.155 | [1.76, 2.36] | 8.11e-40 |
#> | dml plr (rmorie native)  | 2.06***  | 0.15  | [1.76, 2.35] | 2.12e-42 |
#> Consensus (inverse-variance): 2.06 (SE 0.16). holm correction; n = 227.

For a single estimator, morie_estimate_ate() returns one Hajek IPW estimate with its standard error and confidence interval, and morie_estimate_aipw() one augmented IPW estimate:

ipw <- morie_estimate_ate(pair$data, pair$T, pair$Y, pair$covariates)
unlist(ipw[c("ate", "se", "ci_lower", "ci_upper", "n", "ess")])
#>        ate         se   ci_lower   ci_upper          n        ess 
#>   2.052182   0.158635   1.741258   2.363107 227.000000 146.763525

Aggregate companion: incidence-rate ratios

morie_otis_irr_glmm_vm() builds person-year cells from table b01 with the same builder as the a01 design (morie_otis_make_pair_alert_to_volatility_ruhela()) and fits Poisson and negative-binomial GLMs of the regional-volatility count Y_vm_count (region changes) on T_high_ac (alert-state complexity of at least 2), adjusting for year, gender and age band. It has no random effect and no cluster-robust standard error, and says so in its warnings. On the 1,000-row b01 extract, few person-years are treated, so the interval is wide:

irr <- morie_otis_irr_glmm_vm(b01)
irr
#> == IRR Poisson/NB GLM -- vm ~ T_high_ac + demog ==
#>   Cells (id x year): 972
#>   Mean vm / cell: 0.058
#>   T_high_ac prevalence: 1.1%
#>   Treatment T: T_high_ac = 1 if ac >= 2, ac = distinct alert combinations in the person-year (no-alert included)
#>   Outcome Y: vm = Y_vm_count, region changes in the person-year (within-row + row-to-row)
#> 
#> IRR for T_high_ac on vm (covariate-adjusted):
#> Family | IRR | 95% CI | p-value | AIC
#> Poisson | 9.5565 | [4.067, 22.458] | 2.24e-07 | 425.12
#> NegBin2 | 9.5953 | [3.356, 27.436] | 2.46e-05 | 425.01
#> 
#> Interpretation:
#> IRR approx exp(beta) on the alert-complexity treatment. Poisson assumes
#> equidispersion; NB2 relaxes that. Concordance between the two is a
#> robustness signal. All inference is on intra-year (id x FY) cells.
#> 
#> Warnings:
#>   - No random-effect or cluster-robust SE in this implementation; for paper-grade SEs use the OTIS DML pipeline.
#>   - OTIS IDs are year-locked (YYYY-XXXXX-AA); intra-year only.

Mandela classification

Rule 44 of the UN Mandela Rules defines solitary confinement as confinement for 22 hours or more a day without meaningful human contact, and prolonged solitary confinement as solitary confinement for more than 15 consecutive days; Rule 43 prohibits prolonged solitary confinement.

morie_otis_analyze_c11_mandela_classification() applies the 15-day threshold to OTIS table c11 (individuals by aggregate length of confinement). OTIS records no hours out of cell, so the provincial classification is a duration-only proxy:

c11 <- morie_sample("otis_c11")
m11 <- morie_otis_analyze_c11_mandela_classification(c11)
as_table(m11$tables[[1]])
#>   Year                    Type Solitary N Solitary % Torture N Torture % Total
#> 1 2023             Segregation      11071      87.5%      1576     12.5% 12647
#> 2 2023 Restrictive confinement      14231      68.5%      6550     31.5% 20781
#> 3 2024             Segregation       9084      83.5%      1797     16.5% 10881
#> 4 2024 Restrictive confinement      12572      64.0%      7069     36.0% 19641
#> 5 2025             Segregation       7626      79.4%      1982     20.6%  9608
#> 6 2025 Restrictive confinement      14796      59.1%     10249     40.9% 25045

morie_otis_analyze_otis_mandela_provincial_vs_federal() sets those rates beside the federal Structured Intervention Unit figures of Sprott and Doob (2021, Table 19):

pf <- morie_otis_analyze_otis_mandela_provincial_vs_federal(c11)
as_table(pf$tables[[1]])
#>                                Source               Period         Unit
#> 1 Federal SIUs (CSC, Sprott-Doob T19) Nov 2019 - Sept 2020 Person-stays
#> 2      Ontario Provincial Segregation                 2023  Individuals
#> 3      Ontario Provincial Segregation                 2024  Individuals
#> 4      Ontario Provincial Segregation                 2025  Individuals
#> 5     Ontario Provincial RC (broader)                 2023  Individuals
#> 6     Ontario Provincial RC (broader)                 2024  Individuals
#> 7     Ontario Provincial RC (broader)                 2025  Individuals
#>   Solitary % Torture %     N
#> 1      28.4%      9.9%  1960
#> 2      87.5%     12.5% 12647
#> 3      83.5%     16.5% 10881
#> 4      79.4%     20.6%  9608
#> 5      68.5%     31.5% 20781
#> 6      64.0%     36.0% 19641
#> 7      59.1%     40.9% 25045

The units differ (federal person-stays, provincial individuals), as the function’s caveat line states.

Three further classifiers work on single records:

# Federal SIU stay: days, average hours out of cell, % of days missing 4 h
morie_siu_classify_mandela(20, 1.5, 100)$category
#> [1] "Torture"
# OTIS alert flags (mental health, suicide risk, suicide watch) + duration
morie_otis_classify_mandela_combo(1, 1, 0, days = 20)$mandela_category
#> [1] "torture"
# Placement-level rates on b01
mrm_classify_mandela(b01, denominator = "row")
#>     year denominator n_mandela        rate  pct n_broader_rc rate_broader
#> 1   2023         362         0 0.000000000 0.00            0  0.000000000
#> 2   2024         337         3 0.008902077 0.89            3  0.008902077
#> 3   2025         301         4 0.013289037 1.33            4  0.013289037
#> 4 pooled        1000         7 0.007000000 0.70            7  0.007000000

morie_otis_analyze_b05_mandela_classification() applies the same threshold to table b05 (placements by consecutive length). b05 is not in the catalogue, so read it straight from data.ontario.ca (network; not run here):

b05 <- utils::read.csv(paste0(
  "https://data.ontario.ca/dataset/09f7fc65-d3bb-4ca8-8b84-1cdc3ef73c36/",
  "resource/754e8cde-5c74-4c0a-9782-79767a2b26b0/download/",
  "b05_segregation_placements_consecutive_lengths.csv"),
  fileEncoding = "UTF-8-BOM")
m05 <- morie_otis_analyze_b05_mandela_classification(b05)
as_table(m05$tables[[1]])

Federal counterparts

The federal side has its own functions:

names(MORIE_SIUIAP_REPORTS)
#> [1] "final_2024"               "annual_2023_2024"        
#> [3] "preliminary_observations" "thematic_mental_health"  
#> [5] "thematic_indigenous"
morie_siuiap_cite("final_2024")
#> [1] "SIU IAP (2024). SIU IAP Final Report. Public Safety Canada."
doob1 <- analyze_doob_table1_releases()
as_table(doob1$tables[[1]])
#>                Type Revoke (violent) Revoke (non-violent) Revoke (breach)
#> 1        Day Parole      4.8 (0.14%)          33.8 (1.0%)    268.6 (7.9%)
#> 2       Full Parole      5.0 (0.49%)          28.2 (2.8%)     89.4 (8.7%)
#> 3 Statutory Release     84.6 (1.50%)         452.8 (7.8%)  1556.0 (26.7%)
#>       Successful  Total
#> 1 3085.0 (90.9%) 3392.2
#> 2  901.2 (88.0%) 1023.8
#> 3 3735.6 (64.1%) 5829.0

Where to go next