Synchronizing recombinant and parental sequences:

We begin with a master alignment of K putative recombinants and P candidate parental sequences. Represented by a matrix Y = {Y_ns}, , , the alignment is composed of nucleotide base names (A, adenine; G, guanine; T, thymine; C, cytosine) and gap characters (-). To eliminate unnecessary information in Y, we consider only columns where at least one of the recombinants possesses a nucleotide base. For each , we create individual, recombinant-specific alignments Y^(k) by preserving the rows of Y that correspond to recombinant k and its N^(k) − 1 candidate parental sequences (possibly different for each recombinant) and removing the other rows. Sites where the recombinant sequence has a gap are not informative for recombination detection via the DMCP model and are removed from the individual alignments Y^(k). Such gap removal establishes an identity between the lengths of putative recombinants and the number of sites in the individual alignments, , and simplifies information sharing among individual data sets. We map individual alignments onto the master alignment with functions(1)where f_k(i) identifies the site in the master alignment Y that contains the ith nucleotide of recombinant k. Since f_k is a “one-to-one” mapping, f_k(i) ≠ f_k(j) for any i ≠ j, the inverse f_k⁻¹ is defined on the range of f_k that represents the set of sites in the master alignment where the kth recombinant has no deletions.

Dual multiple change-point model:

We assume that conditional on model parameters Φ^(k), each alignment Y^(k) is drawn independently from a DMCP model, i.e.,(2)We first describe the model for evolution of individual alignment sites and then define the model structure across sites.

Columns of each individual alignment Y^(k) are assumed to evolve independently as a continuous-time Markov chain on the state space {A, G, C, T} (Felsenstein 2004). For each site , we parameterize the infinitesimal rate matrix of the Markovian substitution process in terms of its stationary distribution and a transition/transversion rate ratio following Hasegawa et al. (1985). To reduce the number of nuisance parameters in the model, we fix all to the overall observed nucleotide frequencies in Y^(k) (Li et al. 2000). This leaves us with one free parameter defining the substitution matrix . To complete the phylogenetic model specification, we need a bifurcating tree topology describing the historical relationships among nucleotides, with branch lengths representing the expected number of substitutions between the bifurcation events. We further reduce the number of free parameters in the model by integrating out of the likelihood through assuming an exponential prior on each branch length for all . Therefore, the likelihood of site s in recombinant k is a function of three phylogenetic parameters .

To model variation of the phylogenetic parameters along the columns of Y^(k), we assume that the parameters are piecewise constant in s with jumps occurring at unknown change points. We first introduce a set of topology breakpoints , where M^(k) is the unknown number of recombination breakpoints for recombinant k, and , for all . Since topologies can attain only a finite set of values we require that , for all . Similarly we introduce a set of change points for substitution process parameters and assume that and are constant between change points. In summary, our DMCP model for each recombinant k is defined by a set of parameters Φ^(k) = (τ^(k), θ^(k), μ^(k), κ^(k), ρ^(k)), where , , , , and , and the varying dimensionality of the parameter space is determined by M^(k) and J^(k).

Priors for nuisance parameters:

Since our interest in this article is the recombination breakpoints θ^(k), we collect all other parameters for each recombinant into a vector Ψ^(k) = (τ^(k), μ^(k), κ^(k), ρ^(k)) and refer to them as nuisance parameters. We define a prior distribution for nuisance parameters by assuming substantial prior independence, specifically . We assume a noninformative prior for over E^(k) possible tree topologies, relating recombinant k with its potential “parents.” The space of topologies permissible under the DMCP model is formed as described in Minin et al. (2005). Constraints on adjacent topologies are incorporated using a simple Markovian structure(3)The prior distribution for ρ^(k) is specified by first assuming that J^(k) follows a truncated Poisson distribution with a predefined, constant intensity λ and then giving equal prior probabilities to all possible draws of J^(k) integers from the set ,(4)We use one value of λ for all individual alignments as putative recombinant sequences are derived from the same genomic region and therefore should have an approximately equal number of changes in evolutionary pressure. Substitution parameters are a priori log-normally distributed, , , where ν_κ, σ_κ, ν_μ, and σ_μ are either estimated in a hierarchical framework or fixed according to our prior knowledge about sequence variability in the genomic region under study. For more details on specifying the prior distribution for nuisance parameters Ψ^(k) see Minin et al. (2005).

Spatially smoothed prior for recombination locations:

To specify prior probabilities for recombination breakpoint locations, we first switch from their point-process representation to site-specific recombination indicators , where , , , and 1{·} is the indicator function. For clarity of presentation we ignore the fact that the first site of an alignment cannot be a topology breakpoint according to our definition. Such reparameterization allows us to introduce recombination probabilities on the master alignment and then map them onto individual recombinants using functions (1) to define a prior distribution for breakpoint locations,(5)In other words, we determine the prior probability of a site being a recombination breakpoint by finding its position in the master alignment and retrieving the corresponding component from the vector of common recombination probabilities p. Conditional on recombination probabilities p, we assume that breakpoint locations are independent within and between recombinants, so(6)If we denote the number of recombinants that do not have gaps at site s of the master alignment by and define the total number of recombination breakpoints at site s, , for , then Equation 6 simplifies to(7)

Because in practice the total number of observed breakpoints is smaller than the number of sites S by one to two orders of magnitude, estimation of the common recombination probabilities p is unrealistic without further assumptions about their prior distribution. Since HIV recombination is mediated by the enzyme reverse transcriptase that processes nucleotides sequentially (Negroni and Buc 2001), we argue that recombination probabilities should have similar values at adjacent locations. To model such spatial dependency among components of p, we first obtain recombination log odds , where(8)and then use a GMRF prior that penalizes large differences between recombination log odds at neighboring sites,(9)It is easy to see that distribution (9) is improper if we reexpress , where the precision matrix(10)satisfies the identity Q1 = 0. In the context of small area estimation, Ghosh et al. (1998) and Sun et al. (1999) show that despite the singularity of matrix Q such autoregressive priors lead to a proper posterior distribution under mild conditions on the model likelihood function. Our “pseudolikelihood” (7) does not satisfy these conditions when C_s = 0 for all s, or when C_s = T_s for all s. Although very unlikely, such values of recombination counts do have strictly positive probability mass a posteriori. Therefore, we cannot guarantee propriety of the posterior distribution of all model parameters and must replace density (9) by a proper approximation, assuming a priori that , where , I is the S × S identity matrix, and ϵ is a small positive constant. Note that the addition of a positive constant to the diagonal elements of Q preserves the precision matrix sparseness, but forces to be diagonally dominant and therefore positive definite. The proper approximation introduces an additional term, , to the exponent of density (9). In all examples, we use ϵ = 10⁻⁶ such that this term ≈0.05, assuming for all s.

In addition to providing spatial preferences for breakpoint locations, the vector of recombination probabilities p defines the prior distribution for the total number of breakpoints for each alignment k. It is important to put more prior mass on small values of M^(k) to avoid inferring spurious breakpoints from noisy sequence data. The original DMCP model assumes that M^(k) is truncated-Poisson distributed with a rate chosen in such a way that Pr(M^(k) > 0) is equal to a predefined constant, usually 0.5. Similarly, in our hierarchical formulation, we want to control the overall probability of at least one recombination breakpoint in all individual alignments by imposing certain constraints on p. We first note that our site-specific prior on R^(k) imposes a Poisson-binomial distribution for M^(k) with small probabilities of success, usually on the order . Therefore, le Cam's theorem implies that the distribution of M^(k) is approximately Poisson with rate (le Cam 1960). For some constant c, we can set δ_k = −ln(1 − c), so that . Because restricting recombination probabilities for each recombinant individually is impractical, we impose our constraint on the population-level recombination probabilities, . Since , this population-level restriction implies a more conservative prior distribution for the number of breakpoints in each individual data set k with Pr(M^(k) > 0) ≤ c.

We complete our model specification by assuming a priori that ω ∼ Γ(α, β). Following Bernardinelli et al. (1995) we express our prior belief about ω in terms of a ratio of recombination probabilities . On the basis of in vitro HIV recombination detection experiments (Moumen et al. 2001; Dykes et al. 2004; Galetto et al. 2004) we expect that site recombination probabilities should not vary more than sevenfold or equivalently that recombination log odds should not differ by >2. Since our smoothing prior implies that , setting ω = S − 1 ensures that even the most physically distant log odds do not deviate from each other by >2 with probability 0.95. Therefore, we fix the prior mean α/β = S − 1 and choose β to be a small constant (0.01 in the simulation study and 0.02 in the analysis of HIV recombinants).

Inference via MCMC simulation:

To approximate the analytically intractable posterior distribution of all model parameters(11)we sample from (11) using MCMC simulation. During MCMC iterations, we use a Metropolis-within-Gibbs scheme to update the model parameters in two major blocks.

In the first block, we simulate from the full conditional distribution of all individual alignment parameters. The hierarchical structure of our model immediately implies the conditional independence of Φ^(k)s,(12)making it possible to cycle through recombinants for each k and simulate from(13)Minin et al. (2005) describe a reversible-jump MCMC sampler to simulate from the posterior distribution of the DMCP model parameters under a uniform prior on recombination locations. Here, we use a similar algorithm to sample from the distributions in (13) with appropriate modifications of acceptance ratios to incorporate the shared prior over recombination locations. We refer interested readers to Suchard et al. (2003) and Minin et al. (2005) for a more detailed description of the DMCP sampling scheme.

The second block of parameters consists of the recombination log-odds vector γ and the GMRF precision ω. Conditioning on the parameters of the individual alignments yields(14)where recombination counts and trials are as defined in Equation 7, and(15)Note that the sum of recombination probabilities constraint translates into the nonlinear algebraic identity(16)We first describe a sampling procedure on the unconstrained space of recombination log odds and then show how to approximate (16) with a linear constraint that can be incorporated into the sampling algorithm with very little computational burden.

To sample from distribution (14), we rely on the strategy introduced by Rue (2001) and Knorr-Held and Rue (2002) and update (ω, γ) simultaneously. Following their scheme, we first propose a new value for the precision parameter ω* = ωu, where ω is the current precision and u is a random variable with density Pr(u) ∝ 1 + 1/u, defined on the interval [1/U, U], U > 1. This proposal is symmetric and can be tuned by the constant U that controls the “length” of proposal jumps. Given a new value of the precision, we then generate a proposal for the vector of log-odds γ* from a multivariate Gaussian distribution that approximates Pr(γ | C, T, ω*) near its mode, where(17)The Newton–Raphson algorithm is used to locate this mode γ′. Log concavity of density (17) guarantees at most one mode. Then, a second-order Taylor approximation of ln Pr(γ | C, T, ω*) around γ′ generates the proposal mean and precision matrix and concludes the Gaussian proposal construction. The proposed values (ω*, γ*) are accepted or rejected jointly with probability given by the Metropolis–Hastings acceptance ratio. The computational efficiency of this multivariate proposal follows from the special shape of density (17). Note that, during construction of the Gaussian approximation, it is sufficient to apply the Taylor approximation only to the function . Since all mixed derivatives of this function are zero, the off-diagonal elements of are equal to the off-diagonal entries of the Gaussian proposal precision matrix. Therefore, the multivariate normal proposals retain the same sparseness of Q and can be efficiently realized using fast methods of Cholesky decomposition for sparse matrices. For more details on approximating densities of the form similar to (17), see Rue (2001) and Rue et al. (2004).

Implementing prior constraints:

We now turn to the problem of incorporating the imposed restrictions on recombination probabilities into our MCMC algorithm. Implementing a proposal that approximates (17) well while satisfying nonlinear constraint (16) is difficult. However, if we can replace constraint (16) with a linearized form for some vector and scalar e, then we can use unconstrained Gaussian proposals as before to generate a candidate state and recenter the proposal to satisfy the linear constraint via(18)where ν and are obtained via a Taylor expansion of ln Pr(γ | C, T, ω*). Such recentering comes at minimal computational cost as the Cholesky factorization of , needed in the unconditional proposal, can be reused to perform the algebraic operations in Equation 18 (Rue and Held 2005).

To arrive at specific values of a and e, we linearize the function around an arbitrary point . Plugging this linearization into Equation 16 yields(19)(20)Choosing v is less straightforward since we particularly need the linear approximation of (16) to be accurate near a posteriori probable values of γ. To make an intelligent guess about posterior support of recombination log odds, we generate a short “training chain” prior to running our MCMC sampler. During these training iterations, we alternate between sampling from the full conditional distributions defined by (12) and (14) with one heuristic modification that allows us to control the overall recombination probability implicitly. To understand the motivation behind this heuristic, we first point out that , where the latter approximation holds when the number of gaps in the master alignment is small. Therefore, if a new state of γ⁽ⁱ⁾ is accepted at the ith iteration, then a new value of should be close to . To heuristically control the overall recombination probability via the binomial pseudolikelihood (15), we update the vector of trials such that for all s at iteration i, where ⌊x⌋ denotes the largest integer that does not exceed x. After training runs are complete we set the approximation point v to an arithmetic average of simulated components of γ.

In all following analyses we set the prior probability of at least one recombination c = 0.5, and therefore we aim at preserving the condition . Table 1 shows posterior medians and 95% Bayesian credible intervals (BCIs) of the overall recombination probability, , for all analyzed data sets. Although in each case the posterior distribution of the overall recombination is slightly shifted to the right from 0.693, this shift and the spread of the distribution are quite small. Therefore, we conclude that our linear approximation performs well.

Simulation study:

To test our model in the presence of a “known” recombination hotspot, we design a small simulation study. We start with an 888-site long alignment of four primate DNA sequences from humans (H), orangutans (O), squirrel monkeys (S), and lemurs (L), previously used to assess the accuracy of the DMCP model (Minin et al. 2005). This data set strongly supports phylogeny (H, O, (S, L)) as demonstrated by several research groups (Yang and Rannala 1997; Larget and Simon 1999; Suchard et al. 2001). We set the true (under simulation conditions) recombination probabilities for 887 sites (excluding the first site) of this alignment in such a way that sites in the interval [401, 600] are more likely to be breakpoint locations. The sum of recombination probabilities in the hotspot interval is denoted by A. All other sites are assigned a probability of recombination (1 − A)/687, such that recombination probabilities sum to one and define a probability mass function for a discrete random variable attaining values from 2 to 888. We then generate 30 realizations, , of this random variable. For each D_i, , we create a new sequence alignment by permuting the nucleotides of H and L in sites D_i through 888. In these newly formed alignments, sites from 1 to D_i − 1 should support the phylogeny (H, O, (S, L)), while the other portion of the alignment should favor the phylogeny (L, O, (S, H)), obtained by exchanging H and L in the original tree.

We generate only 30 recombinants since this quantity represents well the number of recombinant sequences typically available for analysis. When sample size is small relative to the number of sites covered by putative recombinants, the strength of a hotspot plays a critical role in our ability to recover the region. Therefore, we examine performance of our model for different hotspot probability mass values A = 0.9, 0.7, 0.5, and 0.3. The top left plot in Figure 2 shows the artificially generated probabilities used to simulate recombination events. The remaining plots in the left column of Figure 2 depict posterior medians (solid lines) and 95% BCIs (shaded areas) of recombination probabilities estimated under different true values of hotspot strength A. Solid dots mark true recombination sites where sequences H and L begin their permutation. We see that our model successfully identifies hotspots in the presence of a strong signal. On the other hand, when A = 0.3, simulated breakpoint locations are hardly distinguishable from a random sample from all 888 sites, and our method aptly detects no hotspots. This conservative behavior of our estimation procedure is adequate and, moreover, desirable to avoid erroneous detection of recombination hotspots.

• Download figure
• Open in new tab
• Download powerpoint

Figure 2.—

Simulation study. The left top plot shows recombination probabilities used to simulate recombination events in primate mitochondrial DNA sequences. The letter A denotes the probability mass over the region [401, 600]. The rest of the plots on the left depict the sites at which recombination was simulated (solid dots) and the posterior median (solid line) and 95% BCIs (shading) of inferred site-specific recombination probabilities. The right top plot depicts the prior density of GMRF precision ω. Posterior densities of ω are plotted underneath the prior.

Figure 2, right, shows the prior density (top histogram) and the marginal posteriors of the GMRF precision ω for each value of A. Note that when significant clustering of breakpoints is observed, the posterior mass of ω concentrates closer to zero. We expect such behavior since greater variability in recombination log odds, supported by data, leads to a decrease of smoothness. Additionally, the bottom plot shows that the prior of ω dominates the posterior when true breakpoints are distributed nearly uniformly.

Finally, we test the ability of our model to recover the strength of a hotspot A. For each value of A, Table 2 reports the true proportion of simulated breakpoints contained in the interval [401, 600], the posterior median and 95% BCI of the normalized probability masses, and . The normalization is necessary for comparison of the simulated and estimated recombination probabilities as the former sum to one by construction and the latter sum to ∼ln 2 due to the enforced constraint. Mass consistently underestimates the strength of the hotspot with a 95% BCI covering the true value of A only when A = 0.7. However, if we expand the region by 50 sites in both directions, the posterior of more accurately reflects the true strength of the hotspot. This indicates that uncertainty in estimated breakpoint locations leads to an overestimation of the size of the simulated hotspot region.

A newly observed HIV recombination hotspot:

We apply our model to detect spatial recombination preferences in the gag coding region of the HIV genome. We select 42 sequences from the Los Alamos HIV Sequence Database, all of which have been previously classified as recombinants of pure subtypes A and G (see supplemental information at http://www.genetics.org/supplemental/ for accession numbers). We focus our attention on these two subtypes to limit variation in breakpoint locations due to different subtype composition. Although the effects of such variation remain unknown, experimental evidence is emerging that highlights the importance of subtype composition in the biochemistry of recombination (Chin et al. 2005). The recombinant sequences that we select for our analysis come from several different epidemiological studies (Guo et al. 1993; Durali et al. 1998; Peeters et al. 2000; Barlow et al. 2001; Tebit et al. 2002; Vidal et al. 2003) and therefore should represent a diverse set of recombination events. Besides the recombinant sequences, individual alignments contain representative sequences of subtypes A, G, and B, where the latter serves as an outgroup. Lengths of the alignments range from 562 to 820 nucleotides covering 1118 bp of the HIV genome.

In the top plot of Figure 3 we show the locations of gene products Matrix and Capsid in the master alignment of gag and indicate the position of one of the HIV instability elements (INSs). INSs are RNA sequence motifs involved in post-transcriptional regulation of the HIV gene expression (Mikaélian et al. 1996). It is possible that INS primary or secondary structure promotes recombination. Below the gene map we depict the posterior medians and 95% BCIs of the population-level recombination probabilities. This recombination profile strongly suggests an ∼300-nucleotide-long hotspot near the beginning of the Capsid coding region. The posterior median of the GMRF ω precision amounts to 418, and the 95% BCI of ω is (225, 721).

• Download figure
• Open in new tab
• Download powerpoint

Figure 3.—

Analysis of HIV recombinants. The top plot illustrates the locations of gene products in the HIV gag coding region and marks the position of an instability element (INS) in the Capsid reading frame (hatched box). Below the gene map we show posterior medians (solid line) and 95% BCIs (shading) of population-level recombination probabilities. In the bottom two plots we depict averaged individual-level recombination probabilities (vertical bars), estimated jointly with the hierarchical model (plot second from bottom) and independently with the DMCP model (bottom plot). Solid circles mark breakpoint locations in individual recombinants as estimated by the joint and independent approaches.

The bottom two plots of Figure 3 contain individual-level recombination characteristics, estimated jointly with the hierarchical approach and independently with the DMCP model. In both plots, vertical bars represent naive estimates of site-specific posterior recombination probabilities, obtained by taking the posterior mean of C_s/T_s for all s, where counts C_s and trials T_s retain their definitions in the joint and independent analyses. Solid circles mark point estimates of breakpoint locations in individual alignments. Point estimates are defined as sites where the topology with maximum posterior probability is not equal to the topology with maximum probability at the preceding site. We see that in the joint analysis breakpoints and higher recombination probabilities cluster more tightly in the Capsid region, when compared with the independent DMCP analysis. Moreover, several breakpoints in the “cold” regions of the gag do not receive substantial posterior support during the joint analysis. Under the hierarchical model such shrinkage of individual-level recombination probabilities and breakpoint estimates results from sharing spatial breakpoint information among individual recombinants via the common recombination prior.

A cluster of several breakpoints at the end of Capsid does not substantially elevate the corresponding population-level recombination probabilities. Recombination signal in this region comes only from six recombinants. All of these breakpoints are located at the very end of individual alignments and some of them represent noise, associated with topological uncertainty, rather than recombination events. Figure 3 demonstrates that the posterior support of all breakpoints in this cluster decreases during the joint analysis, resulting in either a shift of their estimates or an elimination of weakly supported breakpoints.

We also compare the joint and independent analyses with respect to their estimates of the total number of breakpoints in each recombinant. We plot the posterior mean number of breakpoints, obtained using both approaches, in Figure 4. Great variability in M^(k) among individual recombinants highlights the importance of allowing flexibility in the number of breakpoints. Most data sets exhibit a slight increase in a posteriori supported number of breakpoints during the joint analysis, but the overall pattern remains unchanged between the two types of analysis. To investigate the cause of the increased support, we compare recombination profiles (data not shown) of all individual recombinants, obtained via the joint and independent approaches. We find that in all cases the increase occurs due to higher values of population-level recombination probabilities in the “hot” portion of gag boosting the posterior confidence in breakpoints located in this region that are weakly or moderately supported under the independent analysis with a flat recombination prior. Therefore, the informative recombination prior does not introduce false breakpoints, but rather amplifies the existing signal inside recombination hotspots. This amplification can be clearly seen by comparing the “skylines” in the bottom two plots of Figure 3.

• Download figure
• Open in new tab
• Download powerpoint

Figure 4.—

Number of breakpoints in individual recombinants. The top plot shows independently estimated posterior mean numbers of breakpoints plotted against jointly estimated posterior mean numbers of breakpoints for the 42 HIV gag individual recombinants. In the two bottom bar plots, we show the posterior mean numbers of breakpoints for the two types of recombination analysis that correspond to the x- and y-axes of the top plot.

Diagnostics of MCMC performance:

To assess the performance of our sampler we first examine parameters at the individual recombinant level of the DMCP models. The total number of breakpoints M^(k) for each alignment k is a pivotal parameter in the DMCP model as its time evolution demonstrates how well our reversible-jump MCMC sampler moves between spaces of different dimension. Since M^(k) is a discrete-valued parameter it is natural to examine the regeneration times t_i, , the time steps at which the Markov chain visits a predefined state (or a set of states), where n is the random number of total visits observed during an MCMC run of fixed length. Mykland et al. (1995) note that the behavior of a renewal process defined by regeneration times of a Markov chain may be used to test the performance of an MCMC sampler. The authors suggest plotting t_i/t_n against i/n. According to the law of large numbers for renewal processes, this scaled regeneration quantile (SRQ) plot should be close to a line passing through points (0, 0) and (1, 1). Since the total number of breakpoints is only a marginalization of the complete Markov chain state, regeneration times of M^(k) are not independent and identically distributed (i.i.d.). However, Li et al. (2000) show that the same interpretation of SRQ plots remains useful even when regeneration times are not strictly i.i.d.

We evaluate the performance of our sampler on the HIV data set with 42 recombinants. For each , we choose the posterior median of M^(k) to be the renewal state for defining regeneration times t_i. We show 42 superimposed SRQ plots in Figure 5, left. Since all SRQ plots in Figure 5 are concentrated around the line y = x, we conclude that our MCMC chains are running long enough to sufficiently sample the posterior distributions of the individual-level parameters.

• Download figure
• Open in new tab
• Download powerpoint

Figure 5.—

Convergence diagnostics. The left plot depicts 42 scaled regeneration quantile (SRQ) plots, where t_i denotes an iteration, at which the total number of breakpoints in individual alignments returns to its posterior median for the ith time. Gelman–Rubin potential scale reduction factors (PSRFs, solid line) and their corresponding 97.5% quantiles (dashed line) for recombination log odds are plotted against site indexes on the right.

To monitor convergence of population-level parameters, we use a Gelman–Rubin potential scale reduction factor (PSRF) (Gelman and Rubin 1992). This statistic tests whether multiple Markov chains, started at different values, converge to the same distribution. The PSRF is approximately equal to the square root of the variance estimated by combining all chains, divided by an average of within-chain variances. If all chains reach stationarity, the PSRF approaches 1. If we assume that the stationary distribution is normal, we also can compute confidence bounds for the t-distributed PSRF.

We calculate PSRFs for the recombination log odds from five chains, each started with different values of γ and ω. We generate initial values ω⁽⁰⁾ from a uniform distribution over the interval (0, 10,000). We then sample a value for from a normal distribution with mean and variance of 2. Conditional on ω⁽⁰⁾ and , we initialize the remaining recombination log odds through a random-walk realization, , for . Such a distribution of starting values for (ω, γ) should be overdispersed with respect to the posterior, as recommended by Gelman and Rubin (1992). Figure 5, right, depicts the PSRFs (solid line) with their 97.5% quantiles for the HIV example recombination log odds. The PSRF and its 95% quantile for ln ω are 1.04078 and 1.09845, respectively. Close proximity of all estimated PSRFs to 1 suggests that all chains reach stationarity.