cosmic

cosmic implements the Conditional Ordinal Stereotype Model for Incident-Level Comparison (COSMIC). It estimates actor-specific propensities from ordinal outcomes observed for multiple actors in the same event. The package was developed to compare police officers' propensities to escalate force, but the model applies more generally when actors share events and the outcome has ordered categories.

COSMIC conditions on the outcomes observed within each incident. This removes incident-level factors shared by the participating officers and supports comparisons of officers relative to peers with whom the data provide informative contrasts. A dynamic programming algorithm evaluates the Poisson-multinomial normalization term, and Stan provides Bayesian inference by Markov chain Monte Carlo.

The methodology is described in:

Ridgeway, G. (2026). A Conditional Ordinal Stereotype Model to Estimate Police Officers' Propensity to Escalate Force. Journal of the American Statistical Association, 1-12

Installation

Install the development version from GitHub with:

install.packages("remotes")
remotes::install_github("gregridgeway/cosmic")

Model fitting uses cmdstanr and CmdStan. If they are not already installed, run:

install.packages(
  "cmdstanr",
  repos = c("https://stan-dev.r-universe.dev", getOption("repos"))
)
cmdstanr::install_cmdstan()

Basic workflow

Input data contain one row for each actor in each event. The event identifier, actor identifier, and ordinal outcome are supplied as columns. Outcomes must be consecutive integers starting at 1.

library(cosmic)

fit <- cosmic(
  data,
  incidentID = incident_id,
  officerID = officer_id,
  y = force_level,
  iter = 2000,
  chains = 4,
  cores = 1,
  threads = 8
)

diagnose(fit)
officers <- officer_summary(fit)
outliers <- outlier_report(officers)

The package provides functions to:

• Fit the COSMIC model with cosmic()
• Check sampler diagnostics with diagnose()
• Extract posterior draws with posterior()
• Compare officers with their data-informed peers using officer_summary()
• Identify officers with strong posterior evidence of being in the upper or lower tail of their peer group using outlier_report()

Demos

The package includes two demos.

Simple six-officer demo

The basic demo simulates 2,000 incidents involving six officers in two disconnected groups. It loads a precomputed fit, checks recovery of the simulated parameters, summarizes officers relative to their peers, and produces an outlier report.

demo("basic_demo", package = "cosmic")

Seattle demo

The Seattle demo analyzes included Seattle Police Department data used in the JASA paper. It fits the full model, checks diagnostics, summarizes officers relative to their peers, and reports officers with strong evidence of being in the tails of their peer distribution. The model fit can take several hours and creates a large fit object.

demo("seattle_demo", package = "cosmic")

For a guided analysis based on precomputed Seattle results, see the package vignette:

vignette("getting-started", package = "cosmic")