2. The Univariate Continuous Phase-Type Distribution: \(T_{MRCA}\) and \(T_{Total}\)
2.1. Theoretical background
In an evolutionary tree the time until two sequences coalesce \(T_i\) can be measured in number of generations \(R_i\) divided by the population size \(N\), this is, \(T_i=R_i/N\). \(T_i\) can easily be proven to approximate to an exponential distribution with rate \(\binom{i}{2}\) (Kingman 1982).
In order to understand the evolutionary history of sequences two additional quantities can be defined –namely the time until the most recent common ancestor \(T_{MRCA}\) and the total tree length \(T_{Total}\). \(T_{MRCA}\) will simply be the sum of all times until two sequences coalesce, in other words \(T_{MRCA}=T_n+T_{n-1}+...+T_2\), where \(T_i\sim\text{exp}(\binom{i}{2})\). \(T_{Total}\), on the other hand, takes into account the length of all possible branches, so \(T_{Total}=nT_n+(n-1)T_{n-1}+...+2T_2\) and, thus, \(iT_i\sim \text{exp}(\frac{i-1}{2})\).
The mean and variance of these two quantities can be derived relatively easily. Defining their distribution, however, is more challenging since both \(T_{MRCA}\) and \(T_{Total}\) are sums of independent exponentially distributed variables with different rates. Their distribution can be computed as a series of convolutions, but their formulation, application and interpretation might be challenging for the average population geneticist.
Instead, we can think of the sum of exponential distributions as a continuous-time Markov chain, where coalescent events are represented as Markov jumps with rate \(T_i\) for \(T_{MRCA}\) and \(iT_i\) for \(T_{Total}\). The Markov chain will end with an absorbing state, which in both cases will be the MRCA.
The Markov chain can be represented using phase-type theory, where the jump rates are defined with a sub-intensity matrix \(T\) and the initial distribution will be defined as a row vector \(\pi\). If we define \(\tau\) as the smallest time (or length) to reach the absorbing state, then \(\tau\sim PH(\pi,T)\). This continuous phase-type distribution has well-documented and easy-to-implement formulas for the expectation, the variance, the survival function, the distribution function and the density function. Moreover, since both \(\pi\) and \(T\) can easily be specified, we can represent evolutionary histories that do not follow the standard coalescent model and still use the same phase-type formulas.
PhaseTypeR contains an efficient implementation of
continuous phase-type distributions. This section shows how to create
phase-type representations of \(T_{MRCA}\) and \(T_{Total}\) under the Kingman’s coalescent
model, and it provides some guidelines for modeling these quantities for
other non-standard coalescent models.
2.2. Time until the most recent common ancestor (\(T_{MRCA}\))
Following Hobolth et al. (2019), we can build a phase-type representation of \(T_{MRCA}\) by defining the sub-intensity matrix, which each row and column will represent the states of the coalescent process.
For example, for Kingman’s coalescent the sub-intensity matrix can be defined as:
\[ \boldsymbol{T} = \begin{bmatrix} -\binom{n}{2} & \binom{n}{2} & 0 & \cdots & 0 \\ 0 & -\binom{n-1}{2} & \binom{n-1}{2} & \cdots & 0 \\ 0 & 0 & -\binom{n-2}{2} & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & 0 & -1 \end{bmatrix} \]
and the initial probability vector
\[ \boldsymbol{\pi} = (1,\ 0, \ ...\,,\ 0) \]
We can therefore define \(T_{MRCA}\)
as following a phase-type distribution using PhaseTypeR
with the following function
# for any integer n > 2
subint_mat_t_mrca <- function(n) {
subint_mat <- diag(-choose(n:2, 2))
subint_mat[-(n-1), -1] <- subint_mat[-(n-1), -1] - subint_mat[-(n-1), -(n-1)]
return(subint_mat)
}For a sample size of \(n=4\), the sub-intensity matrix can be generated using the function above, and the initial probabilities can be defined as a vector of length \(n-1\), starting with one and the rest entries filled with zeros:
n <- 4
subint_mat <- subint_mat_t_mrca(n)
init_probs <- matrix(c(1, rep(0,n-2)), 1, n-1)Now we can use a cont_phase_type class of the
PhaseTypeR package to generate a phase-type representation
of \(T_{MRCA}\) by using
PH():
t_mrca_4 <- PH(subint_mat, init_probs)
t_mrca_4
#> $subint_mat
#> [,1] [,2] [,3]
#> [1,] -6 6 0
#> [2,] 0 -3 3
#> [3,] 0 0 -1
#>
#> $init_probs
#> [,1] [,2] [,3]
#> [1,] 1 0 0
#>
#> $defect
#> [1] 0
#>
#> attr(,"class")
#> [1] "cont_phase_type"There are a number of methods associated with the class
cont_phase_type. For example, the mean for \(T_{MRCA}\) can be computed using:
mean(t_mrca_4)
#> [1] 1.5While the variance can be calculated using:
var(t_mrca_4)
#> [1] 1.1388892.3. Total tree length (\(T_{Total}\)) and the reward-transformation
In a similar way to \(T_{MRCA}\), the total tree length or \(T_{Total}\) can also be represented using phase-type theory (Hobolth, Siri-Jégousse, and Bladt 2019).
For Kingman’s coalescent the sub-intensity matrix can be defined as:
\[ \boldsymbol{T} = \begin{bmatrix} -(n-1)/2 & (n-1)/2 & 0 & \cdots & 0 \\ 0 & -(n-2)/2 & (n-2)/2 & \cdots & 0 \\ 0 & 0 & -(n-3)/2 & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & 0 & -0.5 \end{bmatrix} \]
and the initial probability vector
\[ \boldsymbol{\pi} = (1,\ 0, \ ...\,,\ 0) \]
subint_mat_t_total <- function(n) {
subint_mat = matrix(c(0), nrow = n-1, ncol = n-1)
for (i in 1:n-1) {
subint_mat[i, i] = - 0.5 * (n - i)
if (i < n-1) {
subint_mat[i, i+1] = -subint_mat[i, i]
}
}
subint_mat
}For a sample size of \(n=4\), the sub-intensity matrix and the initial probabilities can be defined as:
n <- 4
subint_mat <- subint_mat_t_total(n)
init_probs <- matrix(c(1, rep(0,n-2)), 1, n-1)\(T_{Total}\) can also be
represented using the cont_phase_type class:
t_total_4 <- PH(subint_mat, init_probs)
t_total_4
#> $subint_mat
#> [,1] [,2] [,3]
#> [1,] -1.5 1.5 0.0
#> [2,] 0.0 -1.0 1.0
#> [3,] 0.0 0.0 -0.5
#>
#> $init_probs
#> [,1] [,2] [,3]
#> [1,] 1 0 0
#>
#> $defect
#> [1] 0
#>
#> attr(,"class")
#> [1] "cont_phase_type"Without further ado, we can use the exact same phase-type formulas
(and thus PhaseTypeR functions) for calculating the mean
and variance of \(T_{Total}\):
mean(t_total_4)
#> [1] 3.666667
var(t_total_4)
#> [1] 5.444444\(T_{MRCA}\) and \(T_{total}\) are tightly interconnected. In fact,
\[ T_{total} = T_{MRCA} \circ [n, n-1, ..., 2] \]
Indeed the core of the sub-intensity matrix for \(T_{MRCA}\) and \(T_{Total}\) is the same. But when \(T_{MRCA}\) measure the time to reach the MRCA, so the length of one branch, \(T_{Total}\) will measure the length of all branches. This leads to proceed to a reward transformation of the original sub-intensity matrix.
The reward-transformation is a way to give more weights to some states than other. So the first state of the sub-intensity matrix should have a rate of \(n\), because before any coalescent event there are \(n\) lineages, then the second state a rate of \(n-1\) and so on.
Therefore, \(T_{total}\) can be represented as a reward-transformed \(T_{MRCA}\), where the reward vector is \([n, n-1, ..., 2]\):
t_total_4_bis <- reward_phase_type(t_mrca_4, c(4, 3, 2))
t_total_4_bis
#> $subint_mat
#> [,1] [,2] [,3]
#> [1,] -1.5 1.5 0.0
#> [2,] 0.0 -1.0 1.0
#> [3,] 0.0 0.0 -0.5
#>
#> $init_probs
#> [,1] [,2] [,3]
#> [1,] 1 0 0
#>
#> $defect
#> [1] 0
#>
#> attr(,"class")
#> [1] "cont_phase_type"Which is exactly the same as our previous definition.
2.4. Density, distribution and quantile functions
PhaseTypeR also includes the density function
(dPH()), quantile function (qPH()),
distribution function (pPH()) and random draw generator
(rPH()) for continuous phase-type distributions. We can
therefore apply these functions for our phase-type representation of
\(T_{MRCA}\) and \(T_{Total}\).
For example, the density function:
dPH(0.5, t_mrca_4)
#> [1] 0.4821092
dPH(1:3, t_total_4)
#> [1] 0.1408529 0.2204939 0.2019975
x <- seq(0, 10, 0.1)
y <- dPH(x, t_mrca_4)
plot(x, y, type = 'l', col = 'orange')
y2 <- dPH(x, t_total_4)
lines(x, y2, col = 'blue')
legend(6, 0.5, legend=c(expression('T'[MRCA]), expression('T'[total])),
col=c("orange", "blue"), lty=1)
title('Density function (n=4)')
The quantile function:
qPH(0.5, t_mrca_4)
#> [1] 1.232831
qPH(c(0.25, 0.75), t_total_4)
#> [1] 1.988284 4.784137
x <- seq(0,0.99,0.01)
y <- qPH(x, t_total_4)
plot(x, y, type = 'l', col = 'blue')
y2 <- qPH(x, t_mrca_4)
lines(x, y2, col = 'orange')
title('Quantile function (n=4)')
legend(0.1, 10, legend=c(expression('T'[MRCA]), expression('T'[total])),
col=c("orange", "blue"), lty=1)
The probability function:
pPH(0.5, t_mrca_4)
#> [1] 0.1214176
pPH(c(0.25, 0.75), t_total_4)
#> [1] 0.001622363 0.030579354
x <- seq(0, 10, 0.1)
y <- pPH(x, t_mrca_4)
plot(x, y, type = 'l', col = 'orange')
y <- pPH(x, t_total_4)
lines(x, y, col = 'blue')
title('Probability function (n=4)')
legend(6, 0.4, legend=c(expression('T'[MRCA]), expression('T'[total])),
col=c("orange", "blue"), lty=1)
And the random drawer:
set.seed(0)
rPH(3, t_mrca_4)
#> [1] 2.1741824 1.5582716 0.9070491
rPH(10, t_total_4)
#> [1] 2.267939 1.919821 1.459578 5.727177 3.080435 3.051581 1.558384 3.619314
#> [9] 1.968178 2.714922set.seed(0)
x <- rPH(10000, t_total_4)
hist(x, main = '10,000 random draws (n=4)', breaks = seq(0, 30, 0.5), ylim = c(0, 3000), xlim=c(0, 20),
col=rgb(0,0,1,0.5))
x <- rPH(10000, t_mrca_4)
hist(x, breaks = seq(0, 30, 0.5), add = T, col=rgb(1,0.5,0,0.5))
legend(15, 2000, legend=c(expression('T'[MRCA]), expression('T'[total])),
col=c("orange", "blue"), lty=1)
box()
