Hamiltonian Monte Carlo with strict convergence criteria reduces run-to-run variability in forensic DNA mixture deconvolution

Taylor D.
Bright J.-A.
Buckleton J.

The interpretation of single source and mixed DNA profiles.

]:

P (V | S_{j}) = \int_{M} P (V | m, S_{j}) f (m) d m

(2)

with the prior probability density function (PDF) of the parameters

f (\cdot)

. In an ideal scenario, LR is independent of any choices made by the laboratory technician and of any random confounding factors. In practice, however, LR depends on the variation in the samples from the crime scene, stochastic events occurring during PCR amplification, allele frequency sampling, parameter settings in the data-processing software,

¹

E.g. GeneMapper™.

hyper-parametrisation of the PG software, etc. [

10

Morrison G.S.

Special issue on measuring and reporting the precision of forensic likelihood ratios: Introduction to the debate.

,

11

Taylor D.
Bright J.-A.
Buckleton J.
Curran J.

An illustration of the effect of various sources of uncertainty on DNA likelihood ratio calculations.

]. Still, even when fixing all of these influences across identical runs on the same EPG data, residual run-to-run variability remains [

8

Perlin M.W.
Legler M.M.
Spencer C.E.
Smith J.L.
Allan W.P.
Belrose J.L.
Duceman B.W.

Validating TrueAllele ® DNA mixture interpretation.

,

Taylor D.
Bright J.-A.
Buckleton J.

The interpretation of single source and mixed DNA profiles.

,

12

Bauer D.W.
Butt N.
Hornyak J.M.
Perlin M.W.

Validating TrueAllele ® interpretation of DNA mixtures containing up to ten unknown contributors.

,

Bright J.-A.
Taylor D.
McGovern C.
Cooper S.
Russell L.
Abarno D.
Buckleton J.

Developmental validation of STRmix™, expert software for the interpretation of forensic DNA profiles.

,

Bright J.-A.
et al.

STRmix™ collaborative exercise on DNA mixture interpretation.

].

So far, this run-to-run variability has been attributed to the inherent stochasticity of the Markov Chain Monte Carlo (MCMC) methods used to estimate mixture model parameters [

Bright J.-A.
et al.

STRmix™ collaborative exercise on DNA mixture interpretation.

,

15

Bright J.-A.
Stevenson K.E.
Curran J.M.
Buckleton J.S.

The variability in likelihood ratios due to different mechanisms.

]. However, as we show here, the apparent run-to-run variability is more likely caused by the choice of convergence criteria used in the MCMC sampler. This is supported by our results demonstrating that run-to-run variability can be reduced when using an MCMC method with strict convergence criteria. For this, we formulate a probabilistic model of DNA mixture deconvolution that only has continuous model parameters, marginalising over the discrete dimensions (genotype sets). While such marginalised models can be properly convergence-controlled, they are generally more expensive to compute. As we show here, however, the intrinsic structure of the posterior probability distribution can be exploited by Hamiltonian Monte Carlo (HMC), maintaining the runtimes of conventional MCMC solutions. We show that the strict convergence criteria afforded by our method significantly reduce run-to-run variability. We further present a deduplication system that efficiently handles the combinatorial growth in the number of genotype sets with increasing numbers of contributors.

1.1 Precision of DNA mixture deconvolution

We define precision of DNA mixture deconvolution as the inverse variance between results of runs with identical hyper-parametrisation on the same EPG data for the same hypotheses. Precision has to be considered in addition to the accuracy of a PG system [

[10]

Morrison G.S.

Special issue on measuring and reporting the precision of forensic likelihood ratios: Introduction to the debate.

], as also the authors of PG algorithms note [

[16]

Gill P.
Benschop C.
Buckleton J.
Bleka Ø.
Taylor D.

A review of probabilistic genotyping systems: EuroForMix, DNAStatistX and STRmix™.

]: “The argument is that the existence of variability [across PG runs — our note] raises doubts about whether any of the results should be accepted.” [

[16]

Gill P.
Benschop C.
Buckleton J.
Bleka Ø.
Taylor D.

A review of probabilistic genotyping systems: EuroForMix, DNAStatistX and STRmix™.

] Courts are often unaware of run-to-run variability, as expert witnesses usually report a single LR number [

[17]

Garrett B.L.
Crozier W.E.
Grady R.

Error rates, likelihood ratios, and jury evaluation of forensic evidence.

]. The issue is even more severe when the verbal scale for reporting LRs is used [

18

Association of Forensic Science Providers
Standards for the formulation of evaluative forensic science expert opinion.

,

19

European Network of Forensic Science Institutes
ENFSI guideline for evaluative reporting in forensic science.

Google Scholar

]. The European Network of Forensic Science Institutes [

[19]

European Network of Forensic Science Institutes
ENFSI guideline for evaluative reporting in forensic science.

Google Scholar

] suggest a scale that defines an LR between 100 and 1000 as “moderately strong support” and LR between 1000 and 10000 as “strong support”. Let us assume that we use a PG system, which, for the given case, outputs results from a normal distribution:

{log}_{10} LR \sim N (2.3, 0.5)

. A single run of the software would give “strong support” in

\approx 73 %

of cases. A technique that reports confidence intervals [

5

Bright J.-A.
Coble M.

Forensic DNA Profiling : A Practical Guide to Assigning Likelihood Ratios.

Google Scholar

,

20

Bright J.-A.
Lee S.-I.
Buckleton J.
Taylor D.

Revisiting the STRmix™ likelihood ratio probability interval coverage considering multiple factors.

Google Scholar

], however, would provide the conservative answer of “moderately strong support”. This highlights the importance of high precision (i.e. low variance) in PG results.

The precision of available commercial solutions has been quantified in several studies [

8

Perlin M.W.
Legler M.M.
Spencer C.E.
Smith J.L.
Allan W.P.
Belrose J.L.
Duceman B.W.

Validating TrueAllele ® DNA mixture interpretation.

,

Taylor D.
Bright J.-A.
Buckleton J.

The interpretation of single source and mixed DNA profiles.

,

12

Bauer D.W.
Butt N.
Hornyak J.M.
Perlin M.W.

Validating TrueAllele ® interpretation of DNA mixtures containing up to ten unknown contributors.

,

Bright J.-A.
Taylor D.
McGovern C.
Cooper S.
Russell L.
Abarno D.
Buckleton J.

Developmental validation of STRmix™, expert software for the interpretation of forensic DNA profiles.

,

Bright J.-A.
et al.

STRmix™ collaborative exercise on DNA mixture interpretation.

]. A standard deviation of LR of

> 1 0^{4}

has been reported between identical runs on a three-contributor mixture (sample ‘3-2’) when the TrueAllele® software was used [

[12]

Bauer D.W.
Butt N.
Hornyak J.M.
Perlin M.W.

Validating TrueAllele ® interpretation of DNA mixtures containing up to ten unknown contributors.

]. Results obtained with the STRmix™ software displayed 10-fold LR difference across runs [

Bright J.-A.
Taylor D.
McGovern C.
Cooper S.
Russell L.
Abarno D.
Buckleton J.

Developmental validation of STRmix™, expert software for the interpretation of forensic DNA profiles.

,

Bright J.-A.
et al.

STRmix™ collaborative exercise on DNA mixture interpretation.

]. In order to increase precision, it has been recommended to increase the number of MCMC iterations, at the expense of a larger computational runtime [

Bright J.-A.
Taylor D.
McGovern C.
Cooper S.
Russell L.
Abarno D.
Buckleton J.

Developmental validation of STRmix™, expert software for the interpretation of forensic DNA profiles.

,

21

Riman S.
Iyer H.
Vallone P.M.

Examining performance and likelihood ratios for two likelihood ratio systems using the PROVEDIt dataset.

].

To determine when to terminate an MCMC sampler in Bayesian inference, convergence criteria are used [

[22]

Gelman A.

Bayesian Data Analysis.

Google Scholar

]. The most popular criterion is the univariate Gelman–Rubin (GR) diagnostic [

22

Gelman A.

Bayesian Data Analysis.

Google Scholar

,

23

Brooks S.P.
Gelman A.

General methods for monitoring convergence of iterative simulations.

Google Scholar

], which compares pooled and within-chain variances of samples to indicate possible convergence. For given model parameters, this diagnostic has a value close to 1 if the samples from different chains result in similar estimates for the marginal distribution. Since actual convergence cannot be quantified as long as the true posterior distribution is not known, convergence criteria measure the stability of samples, and the term “convergence” in MCMC is always relative to the chosen metric.

Some of the available PG software solutions provide users with convergence diagnostics. STRmix™ [

Taylor D.
Bright J.-A.
Buckleton J.

The interpretation of single source and mixed DNA profiles.

] for example calculates the ratio of pooled and within-chain variances of the likelihood of the model (personal communication, Kevin Cheng, Institute of Environmental Science and Research, Ltd., Auckland). GenoProof Mixture [

[24]

Götz F.M.
Schönborn H.
Borsdorf V.
Pflugbeil A.-M.
Labudde D.

Genoproof Mixture 3—new software and process to resolve complex DNA mixtures.

] reports the univariate GRs for the continuous parameters. By default, both software solutions use a predefined constant number of post burn-in samples and then report the value of the diagnostic to the user. If the desired threshold (by default 1.2 in STRmix™, 1.05 in GenoProof Mixture) of the convergence diagnostic has not been achieved, the software offers an option to run additional iterations. The default GR threshold should be rather low. The authors of the diagnostic state [

[22]

Gelman A.

Bayesian Data Analysis.

Google Scholar

]: “The condition of GR near 1 depends on the problem at hand; for most examples, values below 1.1 are acceptable, but for a final analysis in a critical problem, a higher level of precision may be required.” [

[22]

Gelman A.

Bayesian Data Analysis.

Google Scholar

] Providing evidence in court should be considered a critical problem, as the consequences of wrong or doubtful answers are significant [

[25]

Kwong K.

The algorithm says you did it: The use of black box algorithms to analyze complex DNA evidence.

Google Scholar

]. Other scientists researching convergence diagnostics noted [

[26]

Vats D.
Knudson C.

Revisiting the Gelman–Rubin diagnostic.

]: “We argue that a cutoff of GR

\leq

1.1 is much too high to yield reasonable estimates of target quantities.” [

[26]

Vats D.
Knudson C.

Revisiting the Gelman–Rubin diagnostic.

] This seems even more important since the statistical models used in both of these software tools combine continuous (e.g. peak intensities) and discrete (e.g. genotype sets) dimensions. However, GR cannot monitor convergence in discrete dimensions. It is therefore possible that convergence is deduced purely from the continuous parameters, while the genotype set distributions may not have converged at all, offering a possible explanation for the large run-to-run variability observed despite low GR thresholds. In our model, we avoid the issue of assessing convergence of genotype sets by marginalising them out (i.e., the genotype set is not part of the MCMC sampling process).

1.2 Trade-off with execution time

In forensic DNA mixture deconvolution, computational runtime is of great importance, since:

•
Laboratories might have to run software multiple times with different hyper-parametrisations in order to check the robustness of the results or to test hypotheses with different numbers of contributors, different analytical thresholds, different priors, etc.
•
Laboratories might want to quantify the precision of the results over several identical replicates.
•
Forensic laboratories are often working under time pressure, e.g., if cases attract great media attention or laws limit detention time without charges.

In addition to the efficiency of the software implementation, there are multiple factors that influence the execution time, including the number of contributors, the number of alleles per locus, the techniques used to limit the number of considered genotype sets, the convergence criteria, the choice of the optimisation problem (maximum likelihood vs. Bayesian inference), the specification of the model, the priors used, etc.

In general, there is a trade-off between the precision and runtime. Lower runtimes can trivially be achieved by running a smaller number of MCMC iterations, at the expense of precision. Achieving perfect precision is theoretically possible, if runtime is not a priority, by integrating over the latent variables (see Eq. (2)).

The model presented here integrates over the discrete dimensions, i.e., marginalises over genotype sets in order to be able to properly monitor convergence. Such a marginalisation imposes additional computational cost, and a proper estimation of the posterior requires well-suited computational methods. We render our method practically applicable by using a Hamiltonian Monte Carlo (HMC) sampler and by introducing a deduplication system for genotype sets. Thanks to these contributions, we are able to use a strict GR threshold of 1.05 with similar or faster runtime than existing solutions.

2. Materials and methods

We base our probabilistic genotyping model on the work by Taylor et al. [

Taylor D.
Bright J.-A.
Buckleton J.

The interpretation of single source and mixed DNA profiles.

], due to the large number of studies that describe and evaluate this model [

6

Bleka Ø.
Storvik G.
Gill P.

EuroForMix: An open source software based on a continuous model to evaluate STR DNA profiles from a mixture of contributors with artefacts.

,

Taylor D.
Bright J.-A.
Buckleton J.

The interpretation of single source and mixed DNA profiles.

,

27

Taylor D.
Bright J.-A.
Buckleton J.

Using probabilistic theory to develop interpretation guidelines for Y-STR profiles.

,

28

Buckleton J.S.
Lohmueller K.E.
Inman K.
Cheng K.
Curran J.M.
Pugh S.N.
Bright J.-A.
Taylor D.A.

Testing whether stutter and low-level DNA peaks are additive.

,

29

Cheng K.
Bright J.-A.
Kerr Z.
Taylor D.
Ciecko A.
Curran J.
Buckleton J.

Examining the additivity of peak heights in forensic DNA profiles.

Google Scholar

]. The main assumptions in this model are the independence of the loci data given the nuissance parameters and the log-normal distribution of the ratio of the observed EPG peak heights to the peak heights predicted by the generative model (called “expected” peak heights). The generative model consists of several steps, as illustrated in Fig. 1: First, the expected contributions are computed for each genotype set from the considered set of parameters. Then, peak stutter models are applied to predict expected peaks. The next step is to calculate the standard deviation

σ

of the log-normal distribution. Finally, the likelihood of the observed data given the parameters is calculated as the product of the likelihoods of all peaks. The model handles stochastic dropout and drop-in events.

Figure thumbnail gr1 — Fig. 1Illustration of the present probabilistic genotyping model. Top left: the observed peaks $O_{a r}^{l}$ at two selected loci from the green channel of the analysed mixture. Centre left: We analyse the genotype set [(15,16),(14,Q)] for locus D3S1358 and [(6,9),(7,9,3)] for locus TH01. The catch-all dropout allele $Q$ denotes any dropped-out peaks. The lines show the expected contributions for all alleles $t_{a n r}^{l}$ . The function is decreasing within a locus due to degradation modelling. The TH01 locus was modelled with a larger amplification efficiency parameter than D3S1358. Bottom left: Expected peaks $x_{a r}^{l}$ are created by applying stutter ratios to $t_{a n r}^{l}$ for two contributors. In this example, we consider only backward stutter for illustration purposes. Composed peaks (i.e. those that consist of both allelic and stutter contributions) are 14 and 15 at D3S1358, and 6 at TH01. Top right: PDF of the log-normal model for peak 15 at D3S1358 as a function of the expected peak height. Centre right: dropout probability as a function of expected peak height. Bottom right: drop-in probability as a function of observed peak height. Peak 18 at D3S1358 is a drop-in, since it is observed but not expected in the considered genotype set. The dotted lines denote the analytical thresholds.
View Large Image
Figure Viewer
Download Hi-res image
Download (PPT)

We consider the posterior probability of parameters

M

given the evidence

V

P (M | V) = \frac{P (V | M) P (M)}{P (V)} = \frac{P (M) \sum_{j} P (V | M, S_{j}) P (S_{j} | M)}{P (V)} .

(3)

P (S_{j} | M)

is assumed to be a flat prior. This assumption leads to genotype weights that are independent from the person(s) of the interest. In Bayesian inference, evidence is usually neglected as it is too expensive to compute and constant w.r.t. model parameters. We thus obtain the unnormalised posterior

P (M | V) \propto P (M) \sum_{j} P (V | M, S_{j}),

(4)

which we use for estimating

P (V | S_{j})

(see Eq. (2)). We assume peak heights to be conditionally independent given

S_{j}

and

M

, and alleles of a contributor in different loci to be independent from each other. To provide a mathematical formulation of the model, we denote:

•
$O_{a r}^{l}$ : random variable of observed peak height at locus $l$ , allele $a$ , replicate $r$ ;
•
$x_{a r}^{l}$ : expected peak height at locus $l$ , allele $a$ , replicate $r$ ;
•
$f_{X}$ : the PDF of a random variable $X$ ;
•
Q: the “catch-all” dropout allele [
[9]
Taylor D.
Bright J.-A.
Buckleton J.
The interpretation of single source and mixed DNA profiles.
Forensic Sci. Int.: Genetics. 2013; 7: 516-528
Abstract

Full Text

Full Text PDF

PubMed

Scopus (195)

Google Scholar
];
•
$h_{r}^{l}$ : the analytical threshold of the EPG for locus $l$ and replicate $r$ ;
•
$N (0, s^{2})$ : a normal distribution with mean $0$ and standard deviation $s$ ;
•
$c_{ψ}$ : the peak height variance parameter for peak type $ψ$ (e.g. stutter).

The resulting model is a function

x_{a r}^{l} (M, S_{j})

of

M

and

S_{j}

. The likelihood of the observed EPG given parameters

M

and a genotype set

S_{j}

is then:

\sum_{j} P (V | M = m, S_{j}) = \prod_{l} \sum_{j} \prod_{r} P (V_{l, r} | m, S_{j, l}),

(5)

where

V_{l, r}

denotes the available data at locus

l

for replicate

r

, and

S_{j, l}

is the genotype set

j

at locus

l

. The conditional probability is modelled as:

P (V_{l, r} | m, S_{j, l}) = ((\prod_{b \in B} P (O_{b r}^{l} ∣ x_{b r}^{l} (m, S_{j, l}), O_{b r}^{l} > 0)) \prod_{d \in D} P (dropout ∣ x_{d r}^{l} (m, S_{j, l}))) d m .

(6)

The inner multiplications are performed over the set

B

of observed peaks and the set

D

of hypothetical dropout peaks. This model formulates separately the relative likelihood of observed peaks and the probabilities of dropout events. In the following, we abbreviate the notation for

x_{a r}^{l} (m, S_{j, l})

to

x_{a r}^{l}

.

The genotype set weights are calculated by:

P (V | S_{j, l}, l) \propto \sum_{m \in M^{'}} P (m) \prod_{r} P (V_{l, r} | m, S_{j, l}),

(7)

where

M^{'}

is the sample from the posterior.

2.1 Observed peaks

In case a peak is observed (i.e.

O_{a r}^{l} > 0

) our model considers a compound distribution that combines sub-models for peaks that are expected (

f_{Z}

) and for drop-in events

(f_{(O_{a r}^{l} ∣ x_{a r}^{l} = 0, O_{a r}^{l} > 0)})

:

P (O_{a r}^{l} ∣ x_{a r}^{l}, O_{a r}^{l} > h_{r}^{l}) \propto \{\begin{matrix} (1 - d_{rate}) f_{Z} & , if x_{a r}^{l} > 0 \\ d_{rate} f_{(O_{a r}^{l} ∣ x_{a r}^{l} = 0, O_{a r}^{l} > h_{r}^{l})} & , otherwise . \end{matrix}

(8)

For the drop-in events, we use the model introduced by Euroformix [

[6]

Bleka Ø.
Storvik G.
Gill P.

EuroForMix: An open source software based on a continuous model to evaluate STR DNA profiles from a mixture of contributors with artefacts.

]:

x_{a r}^{l} = 0, O_{a r}^{l} > h_{r}^{l} ⟹ O_{a r}^{l} - h_{r}^{l} \sim Exp (λ) .

(9)

Two hyper-parameters based on the level of noise in negative controls are required: the drop-in rate

d_{rate}

and the

λ

of the exponential distribution. For the expected peaks, we assume a log-normal distribution following previous works [

Taylor D.
Bright J.-A.
Buckleton J.

The interpretation of single source and mixed DNA profiles.

]:

Z = ln \frac{O_{a r}^{l}}{x_{a r}^{l}} \sim N (0, σ_{a r}^{l} {(x)}^{2}) .

(10)

The standard deviation of this distribution depends on the components of the expected peaks and their heights:

σ_{a r}^{l} (x) = \frac{1}{x_{a r}^{l}} \sum_{ψ \in Ψ (a)} \frac{c_{ψ} x_{ψ a r}^{l}}{\sqrt{χ_{ψ r}^{l}}} .

(11)

We define

ψ \in Ψ (a) = {a + 2, a + 1, a, a - 1}

. This means that for a single allele

a

, we consider contributions from the allelic peak

a

, the backward stutter from the

a + 1

peak, the forward stutter from the

a - 1

peak, and the double backward stutter from the

a + 2

peak. Our method currently does not consider half-backward stutter.

We use one parameter for allelic peak standard deviation (

c_{ψ} = c_{p}

when

ψ = a

) and a different one for stutter peak standard deviations (

c_{ψ} = c_{s}

when

ψ \neq a

). Individual parameters for the standard deviations of different types of stutter could be introduced without significant changes to the model. The expected peak heights are:

x_{ψ a r}^{l} = \sum_{n} ρ_{(2 a - ψ) a r}^{l} t_{ψ n r}^{l}

x_{a r}^{l} = \sum_{ψ \in Ψ (a)} x_{ψ a r}^{l} .

The product contribution at a selected allele

a

is defined as:

t_{a n r}^{l} = ζ_{a n}^{l} α_{r}^{l} w_{n} e^{0.001 δ_{r n} (m_{a}^{l} - \bar{m})} t_{r}^{'},

(14)

where

m_{a}^{l}

is the molecular weight of the allele, and

\bar{m}

is the average of the largest and smallest observed molecular weights on the EPG within the called peaks. The integer

ζ_{a n}^{l}

is 1 if the genotype of contributor

n

includes allele

a

in locus

l

, and 0 otherwise. It is equal to 2 in case of a homozygote. The scalars

α_{r}^{l}

are the locus-specific amplification efficiencies (LSAE),

w_{n}

are the weights of the contributors that sum up to 1,

δ_{r n}

are the degradation parameters, and

t_{r}^{'}

is the total allelic product expressed in relative fluorescent units (RFU).

Eq. (12) includes the normalised stutter ratios

ρ

and the product contributions from a contributor

n

,

t_{ψ n r}^{l}

. To obtain

ρ

, the unnormalised stutter ratios

π

are deduced from unambiguous profiles (see Supplementary Materials for the list of models with unnormalised stutter ratios). This is done using linear regressions based on allele designation or longest uninterrupted sequences [

https://www.tensorflow.org/probability/api_docs/python/tfp/mcmc/HamiltonianMonteCarlo.

Taylor D.
Bright J.-A.
Buckleton J.

The interpretation of single source and mixed DNA profiles.

]. Normalisation is subsequently required, since multiple types of stutter are considered at the same time:

ρ_{ψ^{'} a r}^{l} = \{\begin{matrix} {(1 + π_{(a - 2) a r}^{l} + π_{(a - 1) a r}^{l} + π_{(a + 1) a r}^{l})}^{- 1} & , if ψ^{'} = a \\ π_{ψ^{'} a r}^{l} ρ_{a a r}^{l} & , otherwise . \end{matrix}

(15)

The sum in Eq. (12) is over the assumed contributors. Finally,

χ_{ψ r}^{l}

models the fact that peak variance is inversely proportional to peak height [

[30]

J.S. Buckleton, The probability of dropout and drop-in, URL https://johnbuckleton.files.wordpress.com/2019/08/the-probability-of-dropout-and-drop-in.pdf.

Google Scholar

]. The rationale behind the formula is explained in Section 1 of the Supplementary Material:

χ_{ψ r}^{l} = \frac{1000}{x_{ψ r}^{l} + 1} + x_{ψ r}^{l} .

(16)

2.2 Dropout events

In case of a dropout, the peak is unobserved because it is below the analytical threshold. This corresponds directly to

P (dropout ∣ x_{a r}^{l}) = P (Z \leq ln \frac{h_{r}^{l}}{x_{a r}^{l}}) .

(17)

See Supplementary Material section 1 for details.

2.3 Parameters of the model and priors

Table 1 summarises the parameters of our model. In typical scenarios, the total number of posterior dimensions will range from 22 (2-contributor, 1-replicate Identifiler™ mixture without modelling of the AMEL locus) to 49 (5-contributor, 3-replicate Fusion 6C mixture with AMEL locus included).

Table 1The parameters of our model. One contributor weight is omitted, as it can be inferred from the other weights (the sum of all weights has to be 1).

L

is the number of loci,

R

is the number of replicates, and

N

is the number of contributors. We further reduce the number of parameters if multiple replicates are performed using the same kit. In such a case, the analysis shares the LSAEs

α_{r}^{l}

across the replicates. The domain boundaries are treated as described in Subsection 3.1 of the Supplementary Materials.

Parameter	Meaning	Domain	Quantity
$c_{s}, c_{p}$	Stutter and peak height standard deviation	$R_{> 0}$	2
$t_{r}^{'}$	Total allelic product	$R_{> 0}$	$R$
$w_{n}$	Weights of contributors	(0,1)	$N - 1$
$δ_{r n}$	Degradation parameters	$R_{> 0}$	$R N$
$α_{r}^{l}$	Locus-specific amplification efficiencies (LSAE)	$R_{> 0}$	$L_{all}$
$i^{2}$	LSAE variance	$R_{> 0}$	1

Open table in a new tab

We define the prior probability of the parameters as:

P (M) = P (α | i) P (i) P (c_{s}) P (c_{p}) .

(18)

Here,

P (α | i)

is the prior ln

(α_{r}^{l}) \sim N (0, i^{2})

that prevents the amplification efficiencies from drifting away from 1 too far. The prior variance

i^{2} \sim Exp (σ_{α})

, where

σ_{α}

is a hyper-parameter to be optimised by the laboratory [

[31]

Bright J.-A.
Taylor D.
Buckleton J.

Statement on minor miscoding in BN formulae STRmix V1.08 with additional comments.

Google Scholar

].

P (c_{s})

and

P (c_{p})

are priors on the peak height standard deviation parameters, which are also present in STRmix™ [

[31]

Bright J.-A.
Taylor D.
Buckleton J.

Statement on minor miscoding in BN formulae STRmix V1.08 with additional comments.

Google Scholar

].

2.4 Other considerations

In order to provide conservative estimates of LR, we use Balding–Nichols sub-population correction [

5

Bright J.-A.
Coble M.

Forensic DNA Profiling : A Practical Guide to Assigning Likelihood Ratios.

Google Scholar

,

32

Balding D.J.
Nichols R.A.

DNA profile match probability calculation: how to allow for population stratification, relatedness, database selection and single bands.

], and we report sub-source LRs [

[33]

Taylor D.
Bright J.-A.
Buckleton J.

The ‘factor of two’ issue in mixed DNA profiles.

] unless specified otherwise. We consider dropout allelic contributions as separate peaks, i.e., (Q,Q) is considered heterozygous. The total number of genotype sets is reduced by considering at most one drop-in per locus and using a drop-in cap. Peaks which are in stutter position and are not included in the genotype set definition are considered drop-ins if they are abnormally tall w.r.t. the origin, see the maximum stutter ratios in Table 3. In summary, if we consider a genotype set for a locus with a tall peak, that genotype set will not be removed during preprocessing if:

•
the corresponding allele score is a part of any of the genotypes from the set;
•
or the height is lower than the drop-in cap and there is no other peak that can be explained as drop-in;
•
or the peak can be explained as a stutter of another peak with an observed stutter ratio lower than the maximum stutter ratio.

2.5 Sampling algorithm

The marginalised (over genotype sets

S_{j}

) PG model from Eq. (4) is prohibitively costly to infer using random-walk MCMC, due to the large number of log-probabilities that need to be calculated when all possible genotype sets are considered. We thus reduce the number of required calculations of log-probabilities by using an adaptive-proposal sampler that has been successfully used in other fields: Hamiltonian Monte Carlo (HMC) [

[34]

Duane S.
Kennedy A.
Pendleton B.J.
Roweth D.

Hybrid Monte Carlo.

]. The difference between HMC and random-walk MCMC is how the proposal distribution is chosen. HMC simulates physical (i.e. Hamiltonian) system dynamics instead of choosing a random point from the neighbourhood of the current sample. This renders HMC very efficient for posteriors with multi-modal or multi-funnel shapes, significant parameter correlations, and/or high dimensionality [

[35]

Betancourt M.

A conceptual introduction to Hamiltonian Monte Carlo.

Google Scholar

]. Unlike GenoProof Mixture, which only changes the value of a single continuous parameter in each iteration, our sampler considers multi-variate moves. The proposal distribution for these moves is dynamically adapted across iterations. It is determined in each iteration based on the local gradient of the log-probability of the model [

[34]

Duane S.
Kennedy A.
Pendleton B.J.
Roweth D.

Hybrid Monte Carlo.

]. This is possible because our model is differentiable, as the discrete genotype sets are marginalised out.

To compute the model gradients, we rely on the proven automatic differentiation framework in TensorFlow Probability [

[36]

Dillon J.V.
Langmore I.
Tran D.
Brevdo E.
Vasudevan S.
Moore D.
Patton B.
Alemi A.
Hoffman M.
Saurous R.A.

TensorFlow distributions.

Google Scholar

], where HMC is also implemented.

²

Thanks to the portability of the TensorFlow library, we can implement both CPU and GPU versions of our inference procedure. However, a naïve TensorFlow implementation would perform below expectations due to the combinatorial growth of the number of possible genotype sets with increasing number of contributors. We therefore introduce an important performance improvement: a deduplication system that can be used with any model that considers all genotypes in a single iteration (e.g. our work, Euroformix). Additional performance optimisations are described in Section 3 of the Supplementary Material.

For the deduplication system, we consider a single locus

l

. The expected peak heights

x_{a r}^{l}

and the standard deviation

σ_{a r}^{l} (x)

of their distribution (Eq. (10)) depend on the genotype sets. If we consider the output for a single allelic position, these values depend only on the continuous parameters and on the alleles at the allelic and stutter (backward, forward, double backward) positions. When multiple genotype sets are considered, the same values are needed multiple times. As an example, consider the locus shown in Fig. 2 and the likely genotype sets from Table 2 for two contributors. Computing every expected peak (and the likelihoods of observing the ratios between the observed and expected peaks) for genotype set {(10,13), (12,15)} entails computations that are also identically required for the other genotype sets. Our deduplication system ensures that each such computation is performed only once, and the result is cached and reused. This leads to large savings in multi-contributor mixtures. As an example, deduplication reduces the number of evaluations of the likelihood ratio of observed and expected peak heights in the D18S51 locus of the MIX13 Case 5 mixture from

\approx

1.3 million to

\approx

86 thousand when 4 unknown contributors are considered in

H_{d}

. We use deduplication during both the calculation of the log-probability and the calculation of the HMC gradient. The deduplication system works as follows: Before a run is performed, we precompute the indices of all duplicate entries in the functions and create two data structures, one containing the information required to evaluate the deduplicated expected peak heights, and one containing the indices for a gather operation to be performed after the log-probabilities of the deduplicated peaks have been calculated. This gather operation then unfolds the deduplicated results into a matrix that stores the log-probabilities per replicate, allelic position, chain, and genotype set. This procedure can be thought of as a cache: the forward propagation first generates the cache by calculating the deduplicated log-probabilities. Then, the cache is used during back-propagation to provide the results for all genotype sets.

Figure thumbnail gr2 — Fig. 2Example of a locus from a DNA mixture with two contributors.
View Large Image
Figure Viewer
Download Hi-res image
Download (PPT)

Table 2Three likely genotype sets

S_{j}

and the resulting contributions to the peaks of the locus in Fig. 2 from two contributors.

C_{n, a} = ρ_{a a 1} t_{a n 1}

is the contribution of contributor

n

from a single copy (i.e.

r = 1

) of allele

a

.

C_{n, a}^{'} = ρ_{(a - 1) a 1} t_{a n 1}

is the contribution from stutter originating from allele

a

. For simplicity, only single-backward stutter and a single replicate are considered. Duplicate entries are highlighted with the same colour.

Table thumbnail fx1001 — Table 2Three likely genotype sets $S_{j}$ **and the resulting contributions to the peaks of the locus in** Fig. 2 **from two contributors.** $C_{n, a} = ρ_{a a 1} t_{a n 1}$ is the contribution of contributor $n$ from a single copy (i.e. $r = 1$ ) of allele $a$ . $C_{n, a}^{'} = ρ_{(a - 1) a 1} t_{a n 1}$ is the contribution from stutter originating from allele $a$ . For simplicity, only single-backward stutter and a single replicate are considered. Duplicate entries are highlighted with the same colour.

3. Results and discussions

Reliable forensic genotyping should minimise the number of false inclusions and false exclusions, exhibit high precision and have low inference runtime. We quantify these aspects for our proposed solution on publicly available test mixtures from previously published benchmarks: the ProvedIt dataset [

[37]

Alfonse L.E.
Garrett A.D.
Lun D.S.
Duffy K.R.
Grgicak C.M.

A large-scale dataset of single and mixed-source short tandem repeat profiles to inform human identification strategies: PROVEDIt.

] and the MIX05 and MIX13 studies [

[38]

Butler J.M.
Kline M.C.
Coble M.D.

NIST interlaboratory studies involving DNA mixtures (MIX05 and MIX13): Variation observed and lessons learned.

]. For the GlobalFiler™ mixtures from the ProvedIt dataset, we use the hyper-parametrisation suggested by Riman et al. [

[21]

Riman S.
Iyer H.
Vallone P.M.

Examining performance and likelihood ratios for two likelihood ratio systems using the PROVEDIt dataset.

] (Table 3). For the MIX05 and MIX13 datasets, we use the hyper-parametrisation from Buckleton et al. [

[39]

Buckleton J.S.
Bright J.-A.
Cheng K.
Budowle B.
Coble M.D.

NIST interlaboratory studies involving DNA mixtures (MIX13): A modern analysis.

]. For MIX05, MIX13, and precision studies we use the FBI extended caucasian population [

[40]

Moretti T.R.
Moreno L.I.
Smerick J.B.
Pignone M.L.
Hizon R.
Buckleton J.S.
Bright J.-A.
Onorato A.J.

Population data on the expanded CODIS core STR loci for eleven populations of significance for forensic DNA analyses in the united states.

] genetic background model. For the challenging mixtures of Section 3.3, we use the NIST 1036-Caucasian background allele frequencies [

[41]

Steffen C.R.
Coble M.D.
Gettings K.B.
Vallone P.M.

Corrigendum to ‘US population data for 29 autosomal STR loci’[Forensic Sci. Int. Genet. 7 (2013) e82–e83].

].

We denote the contributors in the hypotheses by plus-delimited strings. U stands for an unknown contributor, W is a witness, V is a victim, and all other entries denote suspects. The hypothesis V+W+S+U for example has 4 contributors: the victim, the witness, the suspect, and one unknown person. In all benchmark cases, the defendant’s hypothesis is the prosecutor’s hypothesis with the suspect replaced by an unknown contributor. All benchmark mixtures were created in laboratories with known ground-truth genotypes of the contributors.

Table 3Hyper-parameter values used in the present benchmarks. SR stands for stutter ratio. The values of the variance rates and the standard deviation of the locus-specific amplification efficiency (LSAE) were adjusted to our formulation of the model (i.e., natural logarithm instead of

{log}_{10}

, 1/mean for LSAE standard deviation).

Hyper-parameter	ProvedIt	MIX05, MIX13
Drop-in rate ( $d_{rate}$ )	0.0015	0
Drop-in ( $λ$ )	0.032	N/A
Allele variance ( $c_{p}^{2}$ ) shape	5.653	3.57
Allele variance ( $c_{p}^{2}$ ) rate	15.7	5.196
Stutter variance ( $c_{s}^{2}$ ) shape	1.501	6.97
Stutter variance ( $c_{s}^{2}$ ) rate	148.462	9.279
LSAE ( $α_{r}^{l}$ ) standard deviation	32.258	33.333
Rare allele frequency	$2.5 \div 361$	$2.5 \div 202$
Drop-in cap	3 $h_{r}^{l}$
Max. observed backward SR	0.3
Max. observed forward SR	0.15
Max. observed double backward SR	0.05
Wright’s $F_{S T}$ for Balding–Nichols	0.01
Number of chains	4
Number of burnin steps	1200
Leapfrog steps per sample	10
GR stopping threshold	1.05

Open table in a new tab

The linear stutter models are fitted on single-source profiles (forward, backward, and double-backward stutter for ProvedIt and MIX13) or on data provided by the kit manufacturer (only backward stutter for MIX05). The stutter models are available in Section 4 of the Supplementary Material.

All experiments are performed on affordable hardware in the cloud. We use NC8as_T4_v3 instances from Azure Cloud (8 vCPUs, Nvidia Tesla T4 GPU with 16GB VRAM). An exception has been made for ProvedIt Sample 3, which does not fit within 16 GB of VRAM. For this case we rented a Google Cloud virtual machine with a Nvidia A100 GPU with 40GB VRAM.

Table 4Accuracy of our method. Average LRs over 10 runs for our method (“HMC”) in comparison with other solutions on the MIX05 and MIX13 benchmarks with known ground truth. In most cases, all tested methods correctly return LRs larger than 1 for ground-truth inclusion or smaller than 1 for ground-truth exclusion. For MIX13, we compare with the reported results of Euroformix (“EFM”) version 1.11.4 and STRmix™

[39]

Buckleton J.S.
Bright J.-A.
Cheng K.
Budowle B.
Coble M.D.

NIST interlaboratory studies involving DNA mixtures (MIX13): A modern analysis.

. The TPOX locus was ignored in the MIX05 Case 4, due to a tri-allelic pattern of the perpetrator.

Case	$H_{p}$	Truth	${log}_{10} LR$
			HMC	EFM	STRmix™
MIX05 Case 1	Perpetrator+U	Incl.	19.34	–	–
MIX05 Case 2	Perpetrator+U	Incl.	25.67	–	–
MIX05 Case 3	Perpetrator+U	Incl.	22.32	–	–
MIX05 Case 4	Perpetrator+U	Incl.	10.44	–	–
MIX13 Case 1	V+S01A	Incl.	20.15	20.18	20.15
MIX13 Case 2	S02A+U+U	Incl.	17.03	17.28	16.98
MIX13 Case 2	S02B+U+U	Incl.	7.50	7.88	7.26
MIX13 Case 2	S02C+U+U	Incl.	5.41	6.11	5.83
MIX13 Case 2	S02D+U+U	Excl.	−16.18	−2.36	−14.03
MIX13 Case 3	V+W+S03A	Incl.	7.87	6.82	7.69
MIX13 Case 3	V+W+S03B	Excl.	$- \infty$	$- \infty$	$- \infty$
MIX13 Case 4	V+S	Incl.	20.23	19.91	20.15
MIX13 Case 5	S05A+U+U	Incl.	3.38	9.26	3.45
MIX13 Case 5	S05B+U+U	Incl.	1.61	9.38	3.32
MIX13 Case 5	S05C+U+U	Excl.	−8.66	6.45	−9.22
MIX13 Case 5	S05A+U+U+U	Incl.	6.20	–	–
MIX13 Case 5	S05B+U+U+U	Incl.	5.96	–	–
MIX13 Case 5	S05C+U+U+U	Excl.	2.63	–	–

Open table in a new tab

3.1 Accuracy: MIX05 and MIX13 benchmarks

We first benchmark the performance of our method on inter-laboratory studies organised by NIST: MIX05 and MIX13 [

[38]

Butler J.M.
Kline M.C.
Coble M.D.

NIST interlaboratory studies involving DNA mixtures (MIX05 and MIX13): Variation observed and lessons learned.

]. For MIX05, we analyse simultaneously replicates from different kits: ABI’s COFiler, ABI’s SGM Plus, Promega’s Powerplex 16, and ABI’s Profiler Plus. For MIX13, we follow the published studies in using only ABI’s AmpFLSTR IdentiFiler Plus replicate. All cases are analysed with a global analytical threshold of 50 and the ground truth number of contributors with two exceptions: Case 5 from MIX13 is also analysed with 3 contributors (since most laboratories taking part in the original study estimated this number), and Case 2 from MIX13 uses an analytical threshold of 30 (following the recommendation from NIST). For the capillary electrophoresis fragment analysis files, we use default GeneMapper™ ID-X 1.4 analysis settings. The results are presented in Table 4. Our algorithm reproduces most of the results of other solutions, suggesting its validity. Similar to other solutions, our algorithm provides more conservative LR values when a smaller number of contributors is chosen [

[39]

Buckleton J.S.
Bright J.-A.
Cheng K.
Budowle B.
Coble M.D.

NIST interlaboratory studies involving DNA mixtures (MIX13): A modern analysis.

]. The only case in which our model provided LR larger than 1 for a false suspect is S05C in MIX13 Case 5. The genotype of this suspect had been deliberately constructed to share alleles with the true contributors in every locus. 74 out of 108 laboratories have included this suspect in the original study [

[38]

Butler J.M.
Kline M.C.
Coble M.D.

NIST interlaboratory studies involving DNA mixtures (MIX05 and MIX13): Variation observed and lessons learned.

]; our method excludes it in a 3-contributor scenario. For MIX13 Case 4, our algorithm provides a higher LR than the reciprocal of the random match probability. An explanation for this behaviour is given in Subsection 2.1 of the Supplementary Material.

Figure thumbnail gr3 — Fig. 3Precision of our HMC method using different stopping criteria in comparison with STRmix™. We use four different GeneMapper™ analysis methods (A, B, C, D) [14] on Sample 1. Every case is run 100 times, and the resulting per-run LRs are shown as dots. The corresponding maximum likelihood estimations of a normal distribution are shown in the plots below. With standard stopping criteria, our method (cyan) reduces the standard deviation of ${log}_{10} LR$ around 10 fold over STRmix™ (blue). From the published STRmix™ results, we ignored the result provided by participant L4A1 (no known contributor included) and the second run of participant L1A1 (missing in the plots of the original work [14]). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
View Large Image
Figure Viewer
Download Hi-res image
Download (PPT)

3.2 Precision: ProvedIt benchmark

Next we analyse the precision of our method in comparison with STRmix™ on the ProvedIt inter-laboratory benchmark [

Bright J.-A.
et al.

STRmix™ collaborative exercise on DNA mixture interpretation.

]. We use the same analytical thresholds as Kelly et al. [

[42]

Kelly H.
Bright J.-A.
Kruijver M.
Cooper S.
Taylor D.
Duke K.
Strong M.
Beamer V.
Buettner C.
Buckleton J.

A sensitivity analysis to determine the robustness of STRmix™ with respect to laboratory calibration.

] (75, 100, 60, 80, and 100 RFU for the blue, green, yellow, red, and purple dyes, respectively) and report sub-sub-source LRs following Ref. [

Bright J.-A.
et al.

STRmix™ collaborative exercise on DNA mixture interpretation.

].

We focus on the ‘Sample 1’ case, which was previously used to determine the precision of STRmix™ [

Bright J.-A.
et al.

STRmix™ collaborative exercise on DNA mixture interpretation.

].

³

Our method displays high precision also for a simpler ‘Sample 2’ with

{log}_{10} LR = 29.0144 \pm 0.00254

in case of Analysis Method A.

Four different analysis methods of GeneMapper™ are used in comparison (called A, B, C, and D). Analysis method A includes smoothing, no normalisation, and a baseline window of 51. The other methods differ slightly: B uses normalisation, C uses a baseline window of 33, and D does not use smoothing. The results are shown in Fig. 3. For STRmix™, a 10-fold run-to-run variability in the LRs is observed with the default stopping criterion (blue), which has been attributed to the stochastic nature of MCMC [

Bright J.-A.
et al.

STRmix™ collaborative exercise on DNA mixture interpretation.

]. In our HMC method, we check the GR diagnostic every 300 iterations and stop when it is below 1.05 for all parameters. This results in roughly 10 orders of magnitude reduction in the LR standard deviation (cyan). For GeneMapper™ analysis method A, the standard deviation of

{log}_{10} LR

is 10.08 times lower in our method than in STRmix™ [

Bright J.-A.
et al.

STRmix™ collaborative exercise on DNA mixture interpretation.

] (Fig. 3A). For analysis method B (Fig. 3B), it is 8.76 times lower, for analysis method C 10.99 times lower (Fig. 3C), and for analysis method D 10.76 times lower (Fig. 3D).

In the same Fig. 3, we also quantify the influence of the convergence diagnostic on the precision of our algorithm by testing two alternative stopping criteria: In the first, we calculate the mean of GR values and stop when this mean is

<

1.2 (orange). In the second, we use the same convergence metric as STRmix™, but with our threshold value of 1.05 (brown). Due to the efficiency of our HMC sampler, we were unable to simulate chains that resulted in values of this metric reaching as high as 1.2. Interestingly, however, we find that the STRmix™ criterion was sometimes satisfied with threshold 1.05 when some GR values were still

>

1.2.

Taken together, these results show that our approach is able to significantly improve precision over the STRmix™ method, and that this improvement is due to the stricter convergence criteria, as enabled by our purely continuous model parametrisation with genotype sets marginalised out.

Figure thumbnail gr4 — Fig. 4Precision for computationally challenging ProvedIt mixtures. Panel titles indicate the case (top row) and the suspect (bottom row). Each combination is run 10 times, and we plot sub-source ${log}_{10} LR$ . We do not plot exclusion scenarios with LR $=$ 0 for all runs. These are: Sample 2b—non-contributors C18C_Cauc, C99B_AA, and ZT80925_Hisp; Sample 3—non-contributors GT38089_Cauc, JT51484_AA, and WT51359_Cauc. Inclusion cases with correct suspect are plotted in cyan, correct exclusion cases in orange. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
View Large Image
Figure Viewer
Download Hi-res image
Download (PPT)

3.3 Performance: computationally challenging mixtures

Finally, we quantify the performance of our HMC method on the most challenging cases of mixtures with 3 to 5 contributors. This serves to test our approach on cases that are not simple to resolve and that are challenging both from a precision point of view and for computational runtime. For these cases, we use the set of low analytical thresholds from Riman et al. [

[21]

Riman S.
Iyer H.
Vallone P.M.

Examining performance and likelihood ratios for two likelihood ratio systems using the PROVEDIt dataset.

].

We analyse 4 mixtures: Sample 1 and Sample 2 from Bright et al. [

Bright J.-A.
et al.

STRmix™ collaborative exercise on DNA mixture interpretation.

] analysed here without a known contributor (referred here to as Sample 1b and Sample 2b, respectively), and two following 5-person mixtures from ProvedIt with one known contributor:

•
A05_RD14-0003-30_31_32_33_34-1;1;1;1;1-M3I22-0.315GF-Q1.3_01.15 s
(named here Sample 3, known contributor is Contributor 31)
•
E04_RD14-0003-48_49_50_29_30-1;1;2;4;1-M2d-0.279GF-Q2.1_05.15 s
(named here Sample 4, known contributor is Contributor 29)

We use previously suggested analytical thresholds: 35 RFU for the blue channel, 65 for the green, 45 for the yellow, 50 for the red, and 60 for the purple [

[21]

Riman S.
Iyer H.
Vallone P.M.

Examining performance and likelihood ratios for two likelihood ratio systems using the PROVEDIt dataset.

]. The filtered .csv data are used. For each case, we construct all possible prosecutor hypotheses with 1 suspect. We also construct the same number of false hypotheses by choosing the suspects randomly from the NIST 1036 U.S. Population Dataset [

[41]

Steffen C.R.
Coble M.D.
Gettings K.B.
Vallone P.M.

Corrigendum to ‘US population data for 29 autosomal STR loci’[Forensic Sci. Int. Genet. 7 (2013) e82–e83].

]. To quantify precision, we run our method 10 times for each ProvedIt case. The results are shown in Fig. 4. In all cases, our algorithm correctly classifies contributors and non-contributors with high precision. In all but one case where the true contributors are considered, the difference between extremal values of

{log}_{10} LR

is under 0.2. The only exception is Sample 2b when the suspect is Contributor 30. This is further analysed in Section 2 of the Supplementary Material.

The inference runtimes on the benchmark cloud instances are shown in Table 5. We show the results for all the mixtures we analysed with 3 or more unknown contributors. The results are better than the reported runtimes of previous versions of PG software solutions [

[13]

Bright J.-A.
Taylor D.
McGovern C.
Cooper S.
Russell L.
Abarno D.
Buckleton J.

Developmental validation of STRmix™, expert software for the interpretation of forensic DNA profiles.

]. This suggests that despite the increased computational complexity of our marginalised model, the efficiency of HMC sampling and the efficient GPU implementation recover state-of-the-art runtimes as required for practical use of the method.

Table 5Inference runtimes. We report the average (over 10 repetitions) inference times on a single-GPU cloud instance for the listed cases with different numbers of contributors (NoC) in minutes (m) and seconds (s).

Case	NoC	NoC known	Inference time
MIX13 Case 2	3	0	3 m 59 s
MIX13 Case 3	3	0	1 m 59 s
MIX13 Case 5	3	0	4 m 37 s
MIX13 Case 5	4	0	31 m 15 s
ProvedIt Sample 1	4	1	4 m 13 s
ProvedIt Sample 1b	4	0	19 m 39 s
ProvedIt Sample 2b	3	0	6 m 55 s
ProvedIt Sample 3	5	1	59 m 55 s
ProvedIt Sample 4	5	1	27 m 15 s

Open table in a new tab

4. Conclusions and future work

High precision, i.e. low run-to-run variability, of the results provided by probabilistic genotyping methods is key to building trust and to ensuring reliable discriminatory power of the analyses. While run-to-run variability has previously been attributed to the inherent stochasticity of MCMC algorithms [