Fit a hierarchical model of how teams behave

Fits one binomial mixed model per outcome to observed (group, card) decision episodes:

Usage

sdt_fit_behaviour(
  data,
  outcomes = NULL,
  group = "group_id",
  block = "session_id",
  card = "card",
  engine = "lme4",
  nAGQ = 1L,
  frame = c("matched", "all"),
  min_coverage = NULL,
  block_data = NULL,
  block_covariates = NULL,
  verbose = TRUE
)

Arguments

data: A data frame of decision episodes, one row per (group, card).
outcomes: Character vector of outcome names, or NULL to auto-detect.
group, block, card: Column names holding the team identifier, the session (block) the team played in, and the patient card. block may be NULL to fit without a session term. Group labels that repeat across blocks are nested inside their block, with a warning.
engine: "lme4" (default) or "glmmTMB".
nAGQ: Quadrature points for the lme4 engine. 1 (default) is the Laplace approximation. A larger value uses adaptive Gauss-Hermite quadrature, which removes most of the upward bias the Laplace approximation carries in the variance component when clusters hold very few events – the failure mode behind an over-dispersed simulated cohort for a rare behaviour. lme4 supports it only for a model with one grouping factor, so it requires block = NULL; asking for it on a cross-classified model is an error rather than a silent downgrade.
frame: "matched" (default) for the composition-matched frame, or "all" for every row.
min_coverage: Minimum number of groups a card must be seen by to enter the matched frame. NULL chooses it from the data.
block_data: Optional data frame with one row per block, keyed by the block column, holding block-level covariates – typically the output of sdt_course_composition(), which keys on file.
block_covariates: Optional character vector naming the columns of block_data that enter the fixed part of the formula. Passed together with block_data or not at all.
verbose: If TRUE, report progress via cli.

Value

An S3 object of class sdt_behaviour_model:

fits: Named list of fitted models, one per modelled outcome.
outcomes: Names of the outcomes that were modelled.
frame: "matched" or "all", as requested.
terms: The formula template and the column names it was built from, including block_covariates.
data: The frame actually fitted on, with the derived outcome columns attached.
meta: Counts, the frame definition, engine, package versions, runtime, and sds: per outcome, the rows used, the rows dropped as missing, the prevalence, and the group and block standard deviations on the logit scale. With block covariates it also carries composition (the centring and scaling used) and block_variance (the block SD before and after, per outcome).
skipped: One row per outcome not modelled, with the reason.

Details

y ~ card + (1 | block) + (1 | group)

The card fixed effect carries what the casualty is; the group random effect carries what the team is. That second term is the point. It is what lets sdt_simulate_groups() draw a new team rather than another replication of the ones already in the data, and it is what the discrete-event engine in sdt_simulate() has no equivalent of.

The composition-matched frame

Real play is truncated by the clock, so groups see different subsets of the deck. Any between-group statistic computed on that unbalanced frame mixes team behaviour with which cards the team reached, and the second term is much larger than the first. frame = "matched" therefore keeps only cards seen by at least min_coverage groups and only the groups that saw all of them, so every group in the fit answered the same questions. frame = "all" fits on everything and warns: a simulator deals every card to every group, so it cannot reproduce composition spread, and the dispersion ratio is then wrong for a reason that has nothing to do with the model. Measured on the validation corpus, the full deck reports 1.10 for TASC and 1.33 for eFAST where the matched frame built from the same rows reports 1.00 and 1.02.

With min_coverage = NULL the threshold is chosen to maximise cards times groups – the number of complete decision episodes – with ties going to the frame that keeps more groups.

Course composition as a block covariate

block_data is a table with one row per block, keyed by the block column – the shape sdt_course_composition() returns, which keys on file, so fit with block = "file" to join on it. block_covariates names the columns that enter the fixed part, giving

y ~ card + <covariates> + (1 | block) + (1 | group)

The random block intercept stays. The covariates are there to explain part of the session spread, not to replace the term, and with fifteen informative courses a handful of between-course slopes cannot be trusted to replace it.

They also cannot reach further than that. Participants are linked to courses, never to the groups inside them, so composition is constant within a block: it competes with the session random effect and can never explain the team random effect, which is the one sdt_simulate_groups() draws a new team from. A course-level covariate cannot make a simulated team more or less like any other team in its own course.

Numeric and logical covariates are accepted; each is centred and scaled over the distinct blocks, so a course counts once however many groups it ran, and the centring is recorded in $meta$composition$covariates so the simulator can invert it. A covariate that is constant across blocks, that is missing for a block being fitted, or that is collinear with another is refused, and every offender is named in one message.

Each outcome is then fitted twice, with the covariates and without them on exactly the same rows, and $meta$block_variance reports the block SD before and after, the share of block variance removed, and the deviance and AIC change. That doubles the fitting time and it is the point of the exercise: the answer is allowed to be "nothing", and on this corpus it usually is.

Outcomes

outcomes = NULL auto-detects everything the catalogue can build from the columns present: tasc, dcs, etc and conservative from therapy or the did_* flags; efast, xray and ct from did_efast, did_xray and did_ct; red, yellow and green from called_triage; icu and stranded from general_location; ventilated and died from the columns of those names. You may instead pass a character vector naming catalogue outcomes, logical columns of data, or both.

A consensus or modal triage label is deliberately not a source. It is derived from what the teams did, so modelling it would model the model.

What gets skipped

An outcome that cannot be built, is constant, has too few events or groups, or whose fit fails or does not converge is recorded in $skipped with its reason, never dropped in silence. A fit that converges to a zero team variance is kept and flagged singular: that is a finding about the outcome – teams did not differ – and hiding it would hide the negative result the fidelity report exists to surface. A fit whose team SD lands above 5 on the logit scale is flagged boundary for the mirror-image reason: that is separation, and simulated teams will be over-dispersed.

Examples

set.seed(1)
d <- expand.grid(
  card = sprintf("c%02d", seq_len(8)),
  group_id = sprintf("g%02d", seq_len(12)),
  stringsAsFactors = FALSE
)
d$session_id <- ifelse(d$group_id < "g07", "s1", "s2")
b <- stats::rnorm(12, sd = 1.2)
p <- stats::plogis(-0.5 + b[match(d$group_id, unique(d$group_id))])
d$did_efast <- stats::runif(nrow(d)) < p
if (requireNamespace("lme4", quietly = TRUE)) {
  fit <- sdt_fit_behaviour(d, outcomes = "efast", verbose = FALSE)
  fit
}
#> <sdt_behaviour_model> lme4 | matched frame (coverage >= 12)
#> 8 cards x 12 groups in 2 blocks = 96 rows (of 96)
#> 1 outcome modelled in 0.1s:
#>   outcome          prev sd_group sd_block
#>   efast           0.458    1.132    0.000