Note. The Bayesian regularized-horseshoe fit is skipped in the default vignette render to keep R CMD check fast and portable. Set PICMORT_RUN_VIGNETTE_BAYES=true on a host with brms, rstan, and BH installed to include the Bayesian model.

Scope. This vignette mirrors the calibration-first pediatric ICU mortality pipeline end-to-end on a small synthetic cohort (inst/extdata/toy_cohort.rds, n = 80, mortality = 10 %). Numbers shown here are illustrative and do not reproduce the manuscript results. The full pipeline on PIC v1.1.0 (n = 8,736) is shipped as Supplementary File S1 (pipeline/full_pipeline.R) alongside the manuscript bundle.

What this vignette demonstrates

The calibration-first contract for picMort runs through four stages:

Cohort. First ICU stay per patient, age 0–18 y, in-hospital mortality outcome — the frozen contract from vignette("cohort_spec", package = "picMort"). The toy cohort obeys the same schema and the same invariants, downsized to keep the vignette fast.
Features. T₀ + 24 h prediction window, demographics plus vitals and labs aggregated by min / max / mean / last, no length-of-stay, no discharge-time, no post-window leakage. On PIC v1.1.0 this is build_features(cohort, pic_paths()); here, because the toy cohort ships without raw CHARTEVENTS / LABEVENTS, the vignette synthesizes a small, deterministic feature matrix with the same shape and column conventions so the downstream pipeline is exercised faithfully. The synthesis routine is local to this vignette — production runs go through build_features().
Models. Penalized logistic regression (elastic net), XGBoost, and Bayesian logistic regression with a regularized-horseshoe prior. The Bayesian model is the package’s headline functional output: every patient carries a posterior 95 % credible interval on predicted mortality.
Evaluation. Calibration suite (slope, intercept, ICI, calibration-in-the-large) and decision-curve analysis at 5 %, 10 %, 20 % thresholds are the manuscript headline. Discrimination metrics (AUROC, AUPRC, Brier, Brier skill score) are supporting.

library(picMort)
library(data.table)

Stage 1 — load the toy cohort

toy_cohort_path <- system.file("extdata", "toy_cohort.rds", package = "picMort")
cohort <- readRDS(toy_cohort_path)
dim(cohort)
#> [1] 80 13
table(mortality = cohort$hospital_expire_flag)
#> mortality
#>  0  1 
#> 72  8

The toy cohort obeys the same invariants as the production cohort, scaled to a smaller envelope (50–150 rows, mortality 5–20 %).

toy_expectations <- list(
  n_min            = 50L,
  n_max            = 150L,
  mortality_rate   = c(0.05, 0.20),
  age_range_years  = c(0, 18),
  sex_levels       = c("F", "M"),
  distinct_subject = TRUE,
  no_overlap_stays = TRUE
)
assert_cohort_invariants(cohort, expected = toy_expectations)

Stage 2 — build a small illustrative feature matrix

On PIC v1.1.0 the call is features <- build_features(cohort, pic_paths(), window_hours = 24L). The toy cohort ships without raw event tables, so this vignette assembles a matrix of the same shape from the cohort columns plus deterministic synthetic vitals and labs. The hard rule is identical to production: no LOS, no discharge-time, no post-window timestamp ever enters the feature matrix.

make_toy_features <- function(cohort, seed = 20260517L) {
  set.seed(seed)
  n <- nrow(cohort)
  baseline_hr  <- 130 - 4 * pmin(cohort$age_years, 15) + stats::rnorm(n, 0, 8)
  baseline_rr  <- 40  - 1.5 * pmin(cohort$age_years, 15) + stats::rnorm(n, 0, 4)
  baseline_sbp <- 80  + 2.5 * pmin(cohort$age_years, 15) + stats::rnorm(n, 0, 10)
  baseline_spo2 <- pmin(100, 97 + stats::rnorm(n, 0, 1.5))
  baseline_lactate <- pmax(0.5, 1.4 + 0.6 * cohort$hospital_expire_flag +
                             stats::rnorm(n, 0, 0.4))

  x <- data.table::data.table(
    icustay_id          = cohort$icustay_id,
    age_months          = cohort$age_months,
    age_years           = cohort$age_years,
    sex_male            = as.integer(cohort$sex == "M"),
    is_surgical         = as.integer(cohort$is_surgical),
    primary_icd_chapter = cohort$primary_icd_chapter,
    hr_min              = baseline_hr - abs(stats::rnorm(n, 5, 2)),
    hr_max              = baseline_hr + abs(stats::rnorm(n, 8, 3)),
    hr_mean             = baseline_hr,
    hr_last             = baseline_hr + stats::rnorm(n, 0, 4),
    rr_min              = baseline_rr - abs(stats::rnorm(n, 2, 1)),
    rr_max              = baseline_rr + abs(stats::rnorm(n, 3, 1)),
    rr_mean             = baseline_rr,
    rr_last             = baseline_rr + stats::rnorm(n, 0, 2),
    sbp_min             = baseline_sbp - abs(stats::rnorm(n, 6, 2)),
    sbp_max             = baseline_sbp + abs(stats::rnorm(n, 8, 2)),
    sbp_mean            = baseline_sbp,
    sbp_last            = baseline_sbp + stats::rnorm(n, 0, 5),
    spo2_min            = pmin(100, baseline_spo2 - abs(stats::rnorm(n, 2, 1))),
    spo2_max            = pmin(100, baseline_spo2 + abs(stats::rnorm(n, 1, 0.5))),
    spo2_mean           = baseline_spo2,
    spo2_last           = pmin(100, baseline_spo2 + stats::rnorm(n, 0, 1)),
    lactate_min         = baseline_lactate - abs(stats::rnorm(n, 0.2, 0.1)),
    lactate_max         = baseline_lactate + abs(stats::rnorm(n, 0.4, 0.2)),
    lactate_last        = baseline_lactate + stats::rnorm(n, 0, 0.2)
  )

  list(
    x            = x,
    y            = as.integer(cohort$hospital_expire_flag),
    window_hours = 24L,
    feature_set  = "toy"
  )
}

features <- make_toy_features(cohort)
dim(features$x)
#> [1] 80 25

features$x carries the icustay_id join key plus demographics, vitals and a single representative lab. Production runs ship dozens more lab panels; the shape is otherwise identical to build_features() output.

Stage 3 — stratified train / test split

split <- make_train_test_split(features, prop = 0.7, seed = 20260517L)
length(split$train_idx)
#> [1] 55
length(split$test_idx)
#> [1] 25

Stage 4 — fit the three models

Run times below are intentionally aggressive to keep the vignette under two minutes. Production runs use wider hyperparameter grids, 1,000 bootstrap replicates for CIs, and four MCMC chains with 4,000 iterations for the Bayesian fit.

fit_en <- fit_elastic_net(features, split$train_idx, seed = 20260517L)
#> Warning: !  The following column has zero variance so scaling cannot be used:
#>   primary_icd_chapter_symptoms.
#> ℹ Consider using ?step_zv (`?recipes::step_zv()`) to remove those columns
#>   before normalizing.
#> Warning in lognet(x, is.sparse, y, weights, offset, alpha, nobs, nvars, : one
#> multinomial or binomial class has fewer than 8 observations; dangerous ground
#> Warning in lognet(x, is.sparse, y, weights, offset, alpha, nobs, nvars, : one
#> multinomial or binomial class has fewer than 8 observations; dangerous ground
#> Warning in lognet(x, is.sparse, y, weights, offset, alpha, nobs, nvars, : one
#> multinomial or binomial class has fewer than 8 observations; dangerous ground
#> Warning in lognet(x, is.sparse, y, weights, offset, alpha, nobs, nvars, : one
#> multinomial or binomial class has fewer than 8 observations; dangerous ground
#> Warning in lognet(x, is.sparse, y, weights, offset, alpha, nobs, nvars, : one
#> multinomial or binomial class has fewer than 8 observations; dangerous ground
#> Warning in lognet(x, is.sparse, y, weights, offset, alpha, nobs, nvars, : one
#> multinomial or binomial class has fewer than 8 observations; dangerous ground
#> Warning in lognet(x, is.sparse, y, weights, offset, alpha, nobs, nvars, : one
#> multinomial or binomial class has fewer than 8 observations; dangerous ground
#> Warning in lognet(x, is.sparse, y, weights, offset, alpha, nobs, nvars, : one
#> multinomial or binomial class has fewer than 8 observations; dangerous ground
#> Warning in lognet(x, is.sparse, y, weights, offset, alpha, nobs, nvars, : one
#> multinomial or binomial class has fewer than 8 observations; dangerous ground
#> Warning in lognet(x, is.sparse, y, weights, offset, alpha, nobs, nvars, : one
#> multinomial or binomial class has fewer than 8 observations; dangerous ground
#> Warning in lognet(x, is.sparse, y, weights, offset, alpha, nobs, nvars, : one
#> multinomial or binomial class has fewer than 8 observations; dangerous ground
#> Warning in lognet(x, is.sparse, y, weights, offset, alpha, nobs, nvars, : one
#> multinomial or binomial class has fewer than 8 observations; dangerous ground
#> Warning in lognet(x, is.sparse, y, weights, offset, alpha, nobs, nvars, : one
#> multinomial or binomial class has fewer than 8 observations; dangerous ground
#> Warning in lognet(x, is.sparse, y, weights, offset, alpha, nobs, nvars, : one
#> multinomial or binomial class has fewer than 8 observations; dangerous ground
#> Warning in lognet(x, is.sparse, y, weights, offset, alpha, nobs, nvars, : one
#> multinomial or binomial class has fewer than 8 observations; dangerous ground
#> Warning in lognet(x, is.sparse, y, weights, offset, alpha, nobs, nvars, : one
#> multinomial or binomial class has fewer than 8 observations; dangerous ground
#> Warning in lognet(x, is.sparse, y, weights, offset, alpha, nobs, nvars, : one
#> multinomial or binomial class has fewer than 8 observations; dangerous ground
#> Warning in lognet(x, is.sparse, y, weights, offset, alpha, nobs, nvars, : one
#> multinomial or binomial class has fewer than 8 observations; dangerous ground
#> Warning in lognet(x, is.sparse, y, weights, offset, alpha, nobs, nvars, : one
#> multinomial or binomial class has fewer than 8 observations; dangerous ground
#> Warning in lognet(x, is.sparse, y, weights, offset, alpha, nobs, nvars, : one
#> multinomial or binomial class has fewer than 8 observations; dangerous ground
#> Warning in lognet(x, is.sparse, y, weights, offset, alpha, nobs, nvars, : one
#> multinomial or binomial class has fewer than 8 observations; dangerous ground
#> Warning in lognet(x, is.sparse, y, weights, offset, alpha, nobs, nvars, : one
#> multinomial or binomial class has fewer than 8 observations; dangerous ground
#> Warning in lognet(x, is.sparse, y, weights, offset, alpha, nobs, nvars, : one
#> multinomial or binomial class has fewer than 8 observations; dangerous ground
#> Warning in lognet(x, is.sparse, y, weights, offset, alpha, nobs, nvars, : one
#> multinomial or binomial class has fewer than 8 observations; dangerous ground
#> Warning in lognet(x, is.sparse, y, weights, offset, alpha, nobs, nvars, : one
#> multinomial or binomial class has fewer than 8 observations; dangerous ground
#> Warning in lognet(x, is.sparse, y, weights, offset, alpha, nobs, nvars, : one
#> multinomial or binomial class has fewer than 8 observations; dangerous ground
#> Warning in lognet(x, is.sparse, y, weights, offset, alpha, nobs, nvars, : one
#> multinomial or binomial class has fewer than 8 observations; dangerous ground
#> Warning in lognet(x, is.sparse, y, weights, offset, alpha, nobs, nvars, : one
#> multinomial or binomial class has fewer than 8 observations; dangerous ground
#> Warning in lognet(x, is.sparse, y, weights, offset, alpha, nobs, nvars, : one
#> multinomial or binomial class has fewer than 8 observations; dangerous ground
#> Warning in lognet(x, is.sparse, y, weights, offset, alpha, nobs, nvars, : one
#> multinomial or binomial class has fewer than 8 observations; dangerous ground
#> Warning in lognet(x, is.sparse, y, weights, offset, alpha, nobs, nvars, : one
#> multinomial or binomial class has fewer than 8 observations; dangerous ground
#> Warning in lognet(x, is.sparse, y, weights, offset, alpha, nobs, nvars, : one
#> multinomial or binomial class has fewer than 8 observations; dangerous ground
#> Warning in lognet(x, is.sparse, y, weights, offset, alpha, nobs, nvars, : one
#> multinomial or binomial class has fewer than 8 observations; dangerous ground
#> Warning in lognet(x, is.sparse, y, weights, offset, alpha, nobs, nvars, : one
#> multinomial or binomial class has fewer than 8 observations; dangerous ground
#> Warning in lognet(x, is.sparse, y, weights, offset, alpha, nobs, nvars, : one
#> multinomial or binomial class has fewer than 8 observations; dangerous ground
#> Warning in lognet(x, is.sparse, y, weights, offset, alpha, nobs, nvars, : one
#> multinomial or binomial class has fewer than 8 observations; dangerous ground
#> Warning in lognet(x, is.sparse, y, weights, offset, alpha, nobs, nvars, : one
#> multinomial or binomial class has fewer than 8 observations; dangerous ground
c(alpha = fit_en$best_alpha, lambda = fit_en$best_lambda)
#>    alpha   lambda 
#>   0.0000 258.5374

fit_xgb <- fit_xgboost(features, split$train_idx, seed = 20260517L)
#> Warning: !  The following column has zero variance so scaling cannot be used:
#>   primary_icd_chapter_symptoms.
#> ℹ Consider using ?step_zv (`?recipes::step_zv()`) to remove those columns
#>   before normalizing.
fit_xgb$best_nrounds
#> [1] 56

fit_bh <- fit_bayes_horseshoe(
  features, split$train_idx,
  chains = 2L, iter = 500L,
  par_ratio = 0.10, adapt_delta = 0.95,
  seed = 20260517L
)

fit_bh$chains
fit_bh$iter

Stage 5 — predict on the held-out fold

The test fold is touched once. Predictions are aligned by icustay_id; the Bayesian model additionally returns a 95 % credible interval on the predicted probability for every test-fold patient.

x_test <- features$x[split$test_idx, ]
y_test <- features$y[split$test_idx]

pred_en  <- predict_mortality(fit_en,  x_test)
pred_xgb <- predict_mortality(fit_xgb, x_test)
if (!is.null(fit_bh)) {
  pred_bh <- predict_mortality(fit_bh, x_test)
  head(pred_bh, 3)
} else {
  pred_bh <- NULL
}

Stage 6 — evaluation (the manuscript headline)

probs <- list(
  elastic_net = pred_en$prob_raw,
  xgboost     = pred_xgb$prob_raw
)
if (!is.null(pred_bh)) {
  probs$bayes_horseshoe <- pred_bh$prob_raw
}

calib <- lapply(probs, calibration_suite,
                y = y_test, n_boot = 200L, seed = 20260517L)
dca   <- decision_curve(probs, y_test,
                        thresholds = c(0.05, 0.10, 0.20),
                        plot_grid = FALSE)
disc  <- discrimination_metrics(probs, y_test,
                                reference = "elastic_net",
                                n_boot = 200L, seed = 20260517L)

Calibration headline

calib_tbl <- data.table::rbindlist(lapply(names(calib), function(nm) {
  data.table::data.table(
    model     = nm,
    slope     = calib[[nm]]$slope[["estimate"]],
    intercept = calib[[nm]]$intercept[["estimate"]],
    ici       = calib[[nm]]$ici[["estimate"]],
    cit_large = calib[[nm]]$cit_large[["estimate"]]
  )
}))
knitr::kable(calib_tbl, digits = 3,
             caption = "Calibration on the toy-cohort held-out fold.")

Calibration on the toy-cohort held-out fold.
model	slope	intercept	ici	cit_large
elastic_net	701.654	-1.990	0.410	-1.426
xgboost	1.051	-0.324	0.065	-0.234

Decision-curve net benefit at clinically relevant thresholds

dca_headline <- dca[threshold %in% c(0.05, 0.10, 0.20) & type == "model"]
knitr::kable(dca_headline, digits = 4,
             caption = "Net benefit at 5%, 10%, 20% thresholds.")

Net benefit at 5%, 10%, 20% thresholds.
model	threshold	net_benefit	type	lower	upper
elastic_net	0.05	0.0737	model	NA	NA
elastic_net	0.10	0.0222	model	NA	NA
elastic_net	0.20	-0.1000	model	NA	NA
xgboost	0.05	0.0737	model	NA	NA
xgboost	0.10	0.0356	model	NA	NA
xgboost	0.20	0.0400	model	NA	NA

Discrimination (supporting)

knitr::kable(disc, digits = 3,
             caption = "Discrimination metrics with bootstrap 95% CIs.")

Discrimination metrics with bootstrap 95% CIs.
model	metric	estimate	lower	upper
elastic_net	auroc	0.803	0.413	1.000
elastic_net	auprc	0.418	0.057	0.875
elastic_net	brier	0.249	0.249	0.250
xgboost	auroc	0.667	0.043	1.000
xgboost	auprc	0.314	0.022	0.770
xgboost	brier	0.084	0.019	0.157
elastic_net	brier_skill	0.000	0.000	0.000
xgboost	brier_skill	0.663	0.268	0.917

How this maps to the manuscript

This vignette walks the same four-stage pipeline that drives the calibration-first Pediatric Critical Care Medicine manuscript bundle (Supplementary File S1, pipeline/full_pipeline.R). The production bundle:

runs on PIC v1.1.0 (n = 8,736, mortality 0.0844) rather than the 80-row toy cohort,
uses the full feature panel from build_features() (vitals + labs from CHARTEVENTS and LABEVENTS, 13 panels) under the same T₀
- 24 h prediction window,
adds the PIM3 baseline reconstructed via compute_pim3() (Gate G3),
runs the Bayesian fit with four chains and 4,000 iterations, with patient-level 95 % credible intervals as the clinical-functional headline,
reports calibration and decision-curve analysis as the primary outcome, with discrimination supporting.

Numbers in the tables above are not stable across runs of the toy synthesizer and must not be quoted in any paper. The toy cohort exists to let R CMD check exercise the package surface end-to-end and to give new contributors a runnable scaffold.

sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] data.table_1.18.4 picMort_0.1.1     rmarkdown_2.31   
#> 
#> loaded via a namespace (and not attached):
#>   [1] tidyselect_1.2.1      timeDate_4052.112     dplyr_1.2.1          
#>   [4] farver_2.1.2          loo_2.10.0            S7_0.2.2             
#>   [7] fastmap_1.2.0         tensorA_0.36.2.1      rpart_4.1.27         
#>  [10] digest_0.6.39         timechange_0.4.0      lifecycle_1.0.5      
#>  [13] StanHeaders_2.32.10   survival_3.8-6        magrittr_2.0.5       
#>  [16] posterior_1.7.0       compiler_4.6.1        rlang_1.2.0          
#>  [19] sass_0.4.10           tools_4.6.1           yaml_2.3.12          
#>  [22] knitr_1.51            bridgesampling_1.2-1  xgboost_3.2.1.1      
#>  [25] pkgbuild_1.4.8        curl_7.1.0            RColorBrewer_1.1-3   
#>  [28] BH_1.90.0-1           abind_1.4-8           withr_3.0.3          
#>  [31] purrr_1.2.2           sys_3.4.3             nnet_7.3-20          
#>  [34] grid_4.6.1            stats4_4.6.1          sparsevctrs_0.3.6    
#>  [37] future_1.70.0         inline_0.3.21         ggplot2_4.0.3        
#>  [40] MASS_7.3-65           globals_0.19.1        scales_1.4.0         
#>  [43] iterators_1.0.14      cli_3.6.6             mvtnorm_1.4-1        
#>  [46] generics_0.1.4        otel_0.2.0            RcppParallel_5.1.11-2
#>  [49] future.apply_1.20.2   cachem_1.1.0          rstan_2.32.7         
#>  [52] stringr_1.6.0         splines_4.6.1         bayesplot_1.15.0     
#>  [55] parallel_4.6.1        matrixStats_1.5.0     brms_2.23.0          
#>  [58] vctrs_0.7.3           hardhat_1.4.3         V8_8.2.0             
#>  [61] glmnet_5.0            Matrix_1.7-5          jsonlite_2.0.0       
#>  [64] listenv_1.0.0         maketools_1.3.2       foreach_1.5.2        
#>  [67] gower_1.0.2           jquerylib_0.1.4       tidyr_1.3.2          
#>  [70] recipes_1.3.3         glue_1.8.1            parallelly_1.48.0    
#>  [73] codetools_0.2-20      distributional_0.8.1  lubridate_1.9.5      
#>  [76] rsample_1.3.2         stringi_1.8.7         gtable_0.3.6         
#>  [79] shape_1.4.6.1         QuickJSR_1.10.0       tibble_3.3.1         
#>  [82] pillar_1.11.1         furrr_0.4.0           htmltools_0.5.9      
#>  [85] Brobdingnag_1.2-9     ipred_0.9-15          lava_1.9.2           
#>  [88] R6_2.6.1              evaluate_1.0.5        lattice_0.22-9       
#>  [91] backports_1.5.1       bslib_0.11.0          class_7.3-23         
#>  [94] rstantools_2.6.0      Rcpp_1.1.1-1.1        prodlim_2026.03.11   
#>  [97] coda_0.19-4.1         gridExtra_2.3.1       nlme_3.1-169         
#> [100] checkmate_2.3.4       xfun_0.59             buildtools_1.0.0     
#> [103] pkgconfig_2.0.3

Calibration-first pipeline — toy-cohort illustration