Skip to contents

Calibrated statistics with MSU-Fit

Source: vignettes/calibrated-statistics.Rmd
calibrated-statistics.Rmd

A nominal interval is a promise; a calibrated interval is a promise with a measured track record. The MSU-Fit family asks, for every statistical claim the package makes, “when this procedure says 95%, does it cover the truth 95% of the time?” — and measures the answer by Monte Carlo. This guide walks the calibration grid, the verdicts, and the real-data sanity protocol. The estimands being calibrated come from Fitting known signatures.

The three estimands

The calibration suite covers the three statistical claims a signature analysis makes:

  • presence — is this reference signature present in this cohort at all? (the LRT behind ms_test_presence())
  • absolute exposure — given presence, what is the support-conditional interval on the count-scale exposure?
  • zeroing — is a signature truly absent, or just unresolved at this burden? (the compound zeroing claim behind the zero threshold)

For each estimand × generative arm × burden level, the grid simulates catalogs under the known truth, applies the procedure, and records whether the claim held. Coverage that survives the grid is what the verdicts certify.

The calibration grid

ms_calibration_grid() takes the signature dictionary and sweeps a grid of per-sample mutation burdens (n_grid) and true signature shares (shares). m is the number of Monte-Carlo replicates per cell and n_boot the bootstrap size behind the intervals. Both are set small here to keep the guide fast — production calibrations run far more replicates per cell:

ref <- ms_refdb_bundled()$SBS_GRCh38@matrices[[1]]
S   <- ref[, c(1, 5, 9)]
colnames(S) <- paste0("CS", 1:3)
cal <- ms_calibration_grid(S, n_grid = c(200, 1000),
                           shares = c(0.05, 0.2),
                           m = 300, n_boot = 200, seed = 1)
cal
#>      n share cell mean_coverage          se n_reps pass_window true_zero
#> 1  200  0.05    1     0.9200000 0.009043107    300       FALSE         0
#> 2  200  0.20    2     0.9366667 0.008118718    300        TRUE         0
#> 3 1000  0.05    3     0.9422222 0.007777425    300        TRUE         0
#> 4 1000  0.20    4     0.9277778 0.008628517    300       FALSE         0
#>   false_keep false_zero kept_covered kept_missed
#> 1          0        168          679          53
#> 2          0          0          843          57
#> 3          0          1          848          51
#> 4          0          0          835          65

The arm argument switches the generative model between the multinomial and the negative-binomial (overdispersed) arm; a calibration claim only extends to data the arm actually describes.

Verdicts

ms_calibration_verdict() reads the grid and renders the per-estimand verdict against the pre-registered coverage window — whether the procedure may be advertised at its nominal level:

ver <- ms_calibration_verdict(cal)
ver
#> $all_in_window
#> [1] FALSE
#> 
#> $per_cell_in
#> [1] FALSE  TRUE  TRUE FALSE
#> 
#> $failing_cells
#> [1] 1 4
#> 
#> $window
#>   lo   hi 
#> 0.93 0.97

A verdict outside the window is a finding, not a failure of the call: it says the interval family needs recalibration at that operating point, and the package will not paper over it.

The real-data sanity protocol

Calibration on synthetic data raises the obvious question: does the procedure behave on your catalog? ms_calibration_sanity() runs the downsampling sanity on a real catalog: it fits the dictionary at full depth, then at a ladder of dilution levels, and checks that the estimates track the full-depth truth as mutations are removed:

E   <- matrix(c(0.5, 0.3, 0.2), 3, 8,
              dimnames = list(colnames(S), paste0("S", 1:8)))
sim <- ms_simulate(S, exposures = E, burden = 2000, seed = 1)
san <- ms_calibration_sanity(sim, S, n_boot = 100,
                             dilution_sizes = c(500, 1000),
                             n_dilutions = 20, seed = 1)
san
#> $pass
#> [1] FALSE
#> 
#> $s_curve
#>      n         s         se tallies deviation within
#> 1  500 0.8833333 0.01465261     480        NA  FALSE
#> 2 1000 0.9312500 0.01154912     480        NA   TRUE
#> 
#> $synthetic_used
#> [1] FALSE
#> 
#> $n_inversions
#> [1] 0
#> 
#> $inverted_at
#> integer(0)
#> 
#> $systematic_inversion
#> [1] FALSE
#> 
#> $tolerance
#> [1] 0.05
#> 
#> $floor_band
#>  lo  hi 
#> 0.9 1.0 
#> 
#> $share_floor
#> [1] 0.05
#> 
#> $dilution_sizes
#> [1]  500 1000
#> 
#> $n_dilutions
#> [1] 20
#> 
#> $n_boot
#> [1] 100
#> 
#> $bca
#> [1] FALSE
#> 
#> $seed
#> [1] 1
#> 
#> $nb_size
#> [1] 8
#> 
#> $method
#> [1] "likelihood_bidirectional"
#> 
#> $kept
#>    signature sample    truth truth_share
#> 1        CS1     S1 980.7467   0.4903733
#> 2        CS2     S1 602.1586   0.3010793
#> 3        CS3     S1 417.0947   0.2085473
#> 4        CS1     S2 975.5633   0.4877816
#> 5        CS2     S2 619.9607   0.3099803
#> 6        CS3     S2 404.4761   0.2022380
#> 7        CS1     S3 985.3680   0.4926840
#> 8        CS2     S3 567.2798   0.2836399
#> 9        CS3     S3 447.3522   0.2236761
#> 10       CS1     S4 943.2900   0.4716450
#> 11       CS2     S4 619.5991   0.3097995
#> 12       CS3     S4 437.1109   0.2185555
#> 13       CS1     S5 981.8129   0.4909064
#> 14       CS2     S5 608.5195   0.3042598
#> 15       CS3     S5 409.6676   0.2048338
#> 16       CS1     S6 987.1504   0.4935752
#> 17       CS2     S6 610.2330   0.3051165
#> 18       CS3     S6 402.6166   0.2013083
#> 19       CS1     S7 982.3128   0.4911564
#> 20       CS2     S7 557.6302   0.2788151
#> 21       CS3     S7 460.0570   0.2300285
#> 22       CS1     S8 982.9385   0.4914693
#> 23       CS2     S8 608.8342   0.3044171
#> 24       CS3     S8 408.2273   0.2041136
#> 
#> $n_live_samples
#> [1] 8
#> 
#> $dead_samples
#> character(0)
#> 
#> attr(,"class")
#> [1] "ms_calibration_sanity"

The returned table carries the per-signature point estimates and intervals at every dilution level (plus the live/dead sample census) — the “estimates should move smoothly with depth” eyeball check, mechanized. ms_downsample() is the primitive underneath it and is available on its own:

ds <- ms_downsample(sim, depth = 1000, seed = 1)
range(colSums(ds@counts))
#> [1] 1000 1000

Reading the curves

plot_calibration_curve() draws the coverage-vs-burden curves for a calibration grid (shown among the statistical figure types in the visualization gallery). The interpretation rule is the window: inside it the nominal level holds, outside it the operating point is flagged by the verdict.

Where to go next