Getting started with iPEB

Overview

iPEB (improved Parametric Empirical Bayes) analyses longitudinal biomarker data for early-detection screening. It combines a time-gap-aware standardization layer with objective-driven multi-marker weighting: each subject’s history is modelled so that prediction uncertainty grows with the gap between visits, and marker weights are learned to optimize a clinical objective you choose rather than opaque tuning parameters.

library(iPEB)
data(ipeb_example)
str(ipeb_example)
#> 'data.frame':    1439 obs. of  8 variables:
#>  $ id        : int  1 1 1 1 1 1 1 2 2 2 ...
#>  $ case      : int  0 0 0 0 0 0 0 0 0 0 ...
#>  $ time      : num  0 0.865 1.763 2.568 3.406 ...
#>  $ time_to_dx: num  2524 2208 1880 1586 1280 ...
#>  $ m1        : num  8.74 8.4 10.08 9.06 9.57 ...
#>  $ m2        : num  10.9 12 11.1 11.3 11.7 ...
#>  $ m3        : num  8.08 7.91 8.47 8.57 8.78 ...
#>  $ split     : chr  "test" "test" "test" "test" ...

The data are in long format: one row per subject-visit, with a subject id, a case indicator (1 = case, 0 = control), a visit time (years since the first visit), time_to_dx (days from the visit to diagnosis), and biomarkers m1–m3.

train <- subset(ipeb_example, split == "train")
test  <- subset(ipeb_example, split == "test")

Fitting a model

ipeb() fits on the training data. Here we optimize sensitivity at 95% specificity. (We use i.i.d. innovations and no random slope purely to keep this vignette fast; the defaults innovation = "auto" and slope = "auto" choose gap-aware AR(1)/OU innovations and a slope when the data support them.)

fit <- ipeb(train, markers = c("m1", "m2", "m3"),
            objective = "sensitivity", alpha = 0.95,
            innovation = "iid", slope = "off")
fit
#> iPEB model
#>   Objective:          sensitivity
#>   Combiner:           scalar (iPEB-S)
#>   Operating spec (a): 0.95
#>   Detection window:   whole trajectory
#>   Layer:              intercept only, i.i.d. innovations
#>   Markers (3): m1, m2, m3
#>   Weights:
#>     m1  0.919
#>     m2  0.394
#>     m3  -0.023
#>   Training subjects:  168 (56 cases / 112 controls)

The printed summary shows the chosen combiner variant (scalar or multivariate), the operating specificity, the layer configuration, and the learned weights.

Scoring and evaluating new subjects

predict() returns a per-visit iPEB score for new data, and evaluate() reports per-patient sensitivity and median lead time with per-visit specificity at the operating points you request. Thresholds are calibrated on the training controls and applied unchanged to the test subjects.

head(predict(fit, test))
#> [1] -0.06350929  0.11300221  1.34470232 -0.01588068  0.70836508  1.91610740

evaluate(fit, test, specificities = c(0.90, 0.95, 0.99))
#>   specificity sensitivity lead_time realized_specificity       auc
#> 1        0.90   1.0000000 1.2470910            0.8783784 0.9557292
#> 2        0.95   0.9166667 1.1334702            0.9391892 0.9557292
#> 3        0.99   0.8750000 0.3531828            0.9932432 0.9557292

Changing the objective

The same markers can be optimized for a different clinical goal. The lead-time objective rewards earlier detection while retaining sensitivity:

fit_lt <- ipeb(train, markers = c("m1", "m2", "m3"),
               objective = "leadtime", innovation = "iid", slope = "off")
evaluate(fit_lt, test, specificities = 0.95)
#>   specificity sensitivity lead_time realized_specificity     auc
#> 1        0.95   0.9583333  1.193703            0.9425676 0.96875

Feature selection

When a smaller panel is preferred, iPEB can select markers by objective-driven backward elimination to a target size:

fit_sel <- ipeb(train, markers = c("m1", "m2", "m3"),
                objective = "sensitivity", select = "backward", n_markers = 2,
                innovation = "iid", slope = "off")
fit_sel$markers
#> [1] "m1" "m3"

One-call workflow

ipeb_run() fits and evaluates in a single call:

res <- ipeb_run(train, test, markers = c("m1", "m2", "m3"),
                objective = "sensitivity", innovation = "iid", slope = "off",
                specificities = c(0.90, 0.95, 0.99))
res$evaluation
#>   specificity sensitivity lead_time realized_specificity       auc
#> 1        0.90   1.0000000 1.2470910            0.8783784 0.9557292
#> 2        0.95   0.9166667 1.1334702            0.9391892 0.9557292
#> 3        0.99   0.8750000 0.3531828            0.9932432 0.9557292

Notes