knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.align = "center", fig.width = 6, fig.height = 4 ) set.seed(2026)
This vignette walks through the belief-trajectory framework introduced
in mispitools 1.4.0 (Marsico & Egeland, in preparation). The framework
treats the sequential accumulation of per-marker likelihood-ratio evidence
as a trajectory through the probability simplex and computes
diagnostics that are invisible from the combined likelihood ratio alone:
per-step information gain (Bayesian surprise), total-variation path
length, a family of concentration indices, and a per-marker leave-one-out
fragility table.
The useful end-product is the concentration index $C_W^+$, which equals the fraction of the combined weight of evidence carried by the single most impactful supporting marker --- and, by the leave-one-out identity, the fraction of that weight that would be lost under worst-case adversarial removal of that marker. Cases in which $C_W^+$ exceeds the pedigree-specific 90th-percentile simulated under $H_p$ are flagged as fragile inclusions and warrant a leave-one-out review before reporting.
This vignette is organized as an end-to-end walkthrough of a grandparent--grandchild missing-person scenario:
lr_combine().decision_threshold().library(mispitools) library(forrel) library(pedtools)
We use a 2nd-degree grandparent--grandchild pedigree from pedtools.
Under $H_p$ the person of interest is the missing grandchild; under $H_d$
the POI is unrelated to the family.
ped <- linearPed(2) data(Argentina) # Build the allele-frequency list from the Argentina database db <- list() for (m in names(Argentina)[-1]) { freqs <- Argentina[[m]] names(freqs) <- as.character(Argentina$Allele) freqs <- freqs[freqs > 0] db[[m]] <- freqs } marker_names <- names(db)[1:15] db15 <- db[marker_names] for (i in seq_along(db15)) { m <- marker(ped, afreq = db15[[i]], name = names(db15)[i]) if (i == 1) ped <- setMarkers(ped, m) else ped <- addMarkers(ped, m) } ped
We simulate one $H_p$ profile and extract the per-marker likelihood
ratios using forrel::missingPersonLR(). For the vignette we work with a
single profile to keep the exposition compact; in a real calibration
workflow you would simulate thousands and compute the empirical null
distribution of $C_W^+$.
# NOTE: eval=FALSE because the simulation needs forrel's full machinery # which may be slow to run at package-check time. The code below is the # canonical recipe. fullsim <- profileSim(ped, N = 1, ids = labels(ped)) ref <- fullsim for (i in 1:nMarkers(ref)) ref <- setGenotype(ref, id = "3", marker = i, geno = "0/0") poi <- pedtools::singleton("POI") for (i in 1:nMarkers(fullsim)) { af <- afreq(fullsim, marker = i) nm <- name(fullsim, marker = i) mk <- pedtools::marker(poi, afreq = af, name = nm) poi <- if (i == 1) pedtools::setMarkers(poi, mk) else pedtools::addMarkers(poi, mk) } for (i in 1:nMarkers(fullsim)) { g <- pedtools::genotype(fullsim, id = "3", marker = i) poi <- pedtools::setGenotype(poi, id = "POI", marker = i, geno = g) } res <- missingPersonLR(ref, missing = "3", poi = poi, verbose = FALSE) lrs <- setNames(as.numeric(res$LRperMarker), marker_names)
For the remainder of this vignette we use a synthetic per-marker LR
vector that mimics the output of a real missingPersonLR() call. This
lets the vignette run in a few milliseconds at check-time while still
exercising the full trajectory machinery.
# A synthetic matched pair: balanced and concentrated profiles with # similar total log10(LR). lrs_balanced <- c( D3S1358 = 5.2, TH01 = 3.8, D21S11 = 4.1, D18S51 = 2.9, D5S818 = 3.4, D13S317 = 2.7, D7S820 = 4.5, D16S539 = 3.2, CSF1PO = 2.8, vWA = 4.2, TPOX = 3.1, D8S1179 = 3.9, FGA = 5.0, D2S1338 = 3.6, D19S433 = 2.5 ) lrs_concentrated <- c( D3S1358 = 1.3, TH01 = 1.1, D21S11 = 1.4, D18S51 = 1.2, D5S818 = 1.5, D13S317 = 1.2, D7S820 = 1.3, D16S539 = 1.1, CSF1PO = 1.4, vWA = 1.2, TPOX = 1.3, D8S1179 = 1.1, FGA = 250, D2S1338 = 1.5, D19S433 = 1.2 # FGA dominates ) W_balanced <- sum(log10(lrs_balanced)) W_concentrated <- sum(log10(lrs_concentrated)) round(c(W_balanced = W_balanced, W_concentrated = W_concentrated), 2)
Both profiles reach similar total weights of evidence (around 15--17 bans, i.e. a combined LR around $10^{15}$--$10^{17}$). The difference is how the evidence is distributed across markers: the balanced profile has no marker contributing more than about one ban, while the concentrated profile owes almost all its weight to a single strongly-supporting marker (FGA, with $\mathrm{LR} = 250$).
binary_belief_trajectory() turns a per-marker LR vector into a data
frame with one row per step, giving the per-step $\log_{10}\mathrm{LR}$,
the cumulative log-LR, and the posterior $P(H_p)$ at each step.
traj_balanced <- binary_belief_trajectory(lrs_balanced) traj_concentrated <- binary_belief_trajectory(lrs_concentrated) head(traj_balanced)
The full trajectory metrics (entropy, per-step KL divergence, cumulative KL from prior, path length, and the three concentration measures) are computed on the general-form trajectory matrix:
prior <- c(0.5, 0.5) make_traj <- function(lrs) { lr_list <- lapply(unname(lrs), function(r) c(r, 1)) belief_trajectory(prior, lr_list) } metrics_balanced <- trajectory_metrics(make_traj(lrs_balanced)) metrics_concentrated <- trajectory_metrics(make_traj(lrs_concentrated)) summary_df <- data.frame( profile = c("balanced", "concentrated"), path_length = round(c(metrics_balanced$path_length, metrics_concentrated$path_length), 3), concentration = round(c(metrics_balanced$concentration, metrics_concentrated$concentration), 3), herfindahl = round(c(metrics_balanced$concentration_herfindahl, metrics_concentrated$concentration_herfindahl), 3) ) summary_df
The concentrated profile has a concentration index an order of magnitude larger than the balanced one, even though the two reach similar total weights of evidence.
The primary fragility diagnostic is the signed concentration index $C_W^+$, restricted to positive (supporting) contributions. It equals the fractional loss of $W$ under worst-case single-marker removal.
cwp_balanced <- concentration_index_positive(log10(lrs_balanced)) cwp_concentrated <- concentration_index_positive(log10(lrs_concentrated)) round(c(balanced = cwp_balanced, concentrated = cwp_concentrated), 3)
For the concentrated profile, $C_W^+ \approx 0.77$ means the dominant marker (FGA) carries about 77% of the combined weight of evidence. A successful defense challenge of that single marker would leave only about 23% of the combined $W$. The balanced profile has $C_W^+ \approx 0.08$: no single marker carries more than 8% of the total. Both have the same total LR on paper, but they are inferentially very different.
The leave-one-out table makes the per-marker impact explicit:
loo_balanced <- leave_one_out(lrs_balanced) loo_concentrated <- leave_one_out(lrs_concentrated) head(loo_concentrated[order(loo_concentrated$fraction, decreasing = TRUE), ])
The most impactful marker in the concentrated profile is FGA, accounting for the bulk of $W$; all other markers contribute nearly equally small shares. For the balanced profile, no single marker dominates.
The review threshold for $C_W^+$ is not a universal constant; it depends on the laboratory's reference pedigree and marker panel. For the Argentine 15-STR database, the 90th percentile of the $H_p$ null distribution of $C_W^+$ is approximately:
calib <- data.frame( Pedigree = c("linearPed(2) (2nd degree)", "linearPed(3) (3rd degree)", "cousinPed(1) (3rd degree)"), `50%` = c(0.130, 0.140, 0.141), `75%` = c(0.152, 0.163, 0.163), `90%` = c(0.174, 0.192, 0.192), `95%` = c(0.192, 0.211, 0.213), `99%` = c(0.225, 0.254, 0.258), check.names = FALSE ) knitr::kable(calib, caption = "Quantiles of C_W+ under H_p for three pedigrees on the Argentine 15-STR database.")
The concentrated profile above has $C_W^+ \approx 0.77$, which is orders of magnitude above the 99th percentile of any of these pedigree- specific null distributions. Such a case is strongly flagged as fragile: any responsible laboratory workflow should verify the dominant marker (FGA in this example) before reporting, either by re-typing or by independent review. The balanced profile at $C_W^+ \approx 0.08$ is well below even the 50th percentile and passes the fragility filter without flag.
The two steps of the workflow above — calibrating the pedigree-specific
cutoff under $H_p$ and producing a per-case reportable statement — are
exposed as the high-level helpers calibrate_concentration_cutoff() and
fragility_report(). They internally invoke sim_lr_genetic() and
concentration_index_positive(), so a complete pipeline from a
pedtools::ped object to a printable sentence is just two function
calls:
library(pedtools); library(forrel) # Reference pedigree with founder profiles simulated ped <- linearPed(2) ped <- setMarkers(ped, locusAttributes = NorwegianFrequencies[1:15]) ped <- profileSim(ped, N = 1, ids = 2, seed = 1) # Calibrate the pedigree-specific cutoff under H_p (90th percentile by default) cal <- calibrate_concentration_cutoff( reference = ped, missing = 5, numsims = 1500, probs = 0.90, seed = 42 ) cal$cutoff # Per-case fragility report against the calibrated cutoff fr <- fragility_report( per_marker_lrs = lrs_concentrated, cutoff = cal$cutoff, probs = cal$probs ) fr$flag # TRUE -> leave-one-out review required fr$statement # natural-language sentence ready for the case file
The returned statement field is a single character string of the form
"Combined weight of evidence W = w bans, C_W+ = c, ... a successful
challenge of the top marker would leave (1-c)w bans of residual support",
suitable to paste verbatim into the laboratory's case report.
Real missing-person casework often combines DNA evidence with
non-genetic evidence (sex, age, phenotypic traits). mispitools
implements the Egeland--Marsico (2026) Markov-chain framework for
computing per-item supplementary-evidence LRs, and extends the trajectory
framework to combined (FDE + SE) sequences.
# Supplementary evidence LRs (numerical example from Egeland & Marsico 2026) LR_sex <- 2.0 # lr_sex(...) LR_age <- 7.7 # lr_age(...) LR_hair <- 3.3 # lr_hair_color(...) # Extended per-marker LR vector: 15 STRs + 3 SE items lrs_combined <- c(lrs_balanced, Sex = LR_sex, Age = LR_age, Hair = LR_hair) traj_combined <- binary_belief_trajectory(lrs_combined) # Trajectory metrics and concentration on the combined 18-step sequence metrics_combined <- trajectory_metrics(make_traj(lrs_combined)) cwp_combined <- concentration_index_positive(log10(lrs_combined))
The combined trajectory has the same structure as a pure-STR trajectory: each SE item is simply an additional step. The concentration index computed on the combined 18-step sequence reflects the most impactful contributor overall, whether it is a STR or an SE item.
Once you have per-marker LRs for a case, you can feed the combined $W$
into mispitools's decision-threshold machinery (decision_threshold(), plot_decision_curve(),
LRdist()) to obtain a weighted-error-minimized cutoff. The fragility
diagnostic of $C_W^+$ complements the decision threshold: a case may
exceed the inclusion threshold while still being flagged as fragile, in
which case the analyst should perform a leave-one-out review before
reporting the inclusion.
?belief_trajectory,
?trajectory_metrics, ?concentration_index_positive,
?leave_one_out, and ?familias_trajectory for the full API.Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.