Processing math: 81%

Tracing the likelihood calculation of a Gaussian model

Venelin Mitov

2024-03-15

This tutorial shows how the pruning likelihood calculation algorithm works for GLInv models. This is an advanced topic, which can be helpful in case you wish to validate another implementation of the algorithm, e.g. in a different programming language, such as Java. For details about the mathematical terms and equations involved in the following calculations, refer to (Mitov et al. 2019).

Example tree and data

library(ape); 
library(PCMBase);

# Non-ultrametric phylogenetic tree of 5 tips in both examples:
treeNewick <- "((5:0.8,4:1.8)7:1.5,(((3:0.8,2:1.6)6:0.7)8:0.6,1:2.6)9:0.9)0;"
tree <- PCMTree(read.tree(text = treeNewick))
# Partitioning the tree in two parts and assign the regimes:
PCMTreeSetPartRegimes(tree, part.regime = c(`6`=2), setPartition = TRUE, inplace = TRUE)

pOrder <- c(PCMTreeGetLabels(tree)[tree$edge[PCMTreePostorder(tree), 2]], "0")

# Trait-data:
X <- cbind(
  c(0.3, NaN, 1.4), 
  c(0.1, NaN, NA), 
  c(0.2, NaN, 1.2), 
  c(NA, 0.2, 0.2), 
  c(NA, 1.2, 0.4))

colnames(X) <- as.character(1:5)

A tree with five tips and two evolutionary regimes

A mixed Gaussian model

model.OU.BM <- MixedGaussian(
  k = nrow(X), 
  modelTypes = c(
    BM = "BM__Omitted_X0__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Sigma_x__Omitted_Sigmae_x",
    OU = "OU__Omitted_X0__H__Theta__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Sigma_x__Omitted_Sigmae_x"), 
  mapping = c(2, 1), 
  Sigmae_x = structure(
    0, 
    class = c("MatrixParameter", "_Omitted", 
              description = "upper triangular factor of the non-phylogenetic variance-covariance")))

model.OU.BM <- PCMApplyTransformation(model.OU.BM)
model.OU.BM$X0[] <- c(NA, NA, NA)
model.OU.BM$`1`$H[,,1] <- cbind(
  c(.1, -.7, .6), 
  c(1.3, 2.2, -1.4), 
  c(0.8, 0.2, 0.9))
model.OU.BM$`1`$Theta[] <- c(1.3, -.5, .2)
model.OU.BM$`1`$Sigma_x[,,1] <- cbind(
  c(1, 0, 0), 
  c(1.0, 0.5, 0), 
  c(0.3, -.8, 1))

model.OU.BM$`2`$Sigma_x[,,1] <- cbind(
  c(0.8, 0, 0), 
  c(1, 0.3, 0), 
  c(0.4, 0.5, 0.3))

print(
  PCMTable(model.OU.BM, removeUntransformed = FALSE), 
  xtable = TRUE, type=tableOutputType)


regime type X0 H Θ Σx Σ
:global: NA [...]
1 OU [0.101.300.80−0.702.200.200.60−1.400.90] [1.30−0.500.20] [1.001.000.300.000.50−0.800.000.001.00] [2.090.260.300.260.89−0.800.30−0.801.00]
2 BM [0.801.000.400.000.300.500.000.000.30] [1.800.500.120.500.340.150.120.150.09]

Likelihood values for the three example variants

options(digits = 4)
# Variant 1: 
PCMLik(X[, tree$tip.label], tree, model.OU.BM)
## [1] -11.92
## attr(,"X0")
## [1]  9.566 -6.349 15.254
# Variant 2: First we call the function PCMInfo to obtain a meta-information object.
metaI.variant2 <- PCMInfo(X[, tree$tip.label], tree, model.OU.BM)
# Then, we manually change the vector of present coordinates for the root node.
# The pc-matrix is a k x M matrix of logical values, each column corresponding
# to a node. The active coordinates are indicated by the TRUE entries.
# To prevent assigning to the wrong column in the pc-table, we first assign
# the node-labels as column nanmes.
colnames(metaI.variant2$pc) <- PCMTreeGetLabels(tree)
metaI.variant2$pc[, "0"] <- c(TRUE, FALSE, TRUE)
# After the change, the pc-matrix looks like this:
metaI.variant2$pc
##          5     4     3     2     1     0    7     9     8     6
## [1,] FALSE FALSE  TRUE  TRUE  TRUE  TRUE TRUE  TRUE  TRUE  TRUE
## [2,]  TRUE  TRUE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE
## [3,]  TRUE  TRUE  TRUE FALSE  TRUE  TRUE TRUE  TRUE  TRUE  TRUE
# And the log-likelihood value is:
PCMLik(X[, tree$tip.label], tree, model.OU.BM, metaI = metaI.variant2)
## [1] -12.17
## attr(,"X0")
## [1] 10.611    NaN  8.747
# Variant 3: We set all NaN values in X to NA, to indicate that these are
# missing measurements
X3 <- X
X3[is.nan(X3)] <- NA_real_
PCMLik(X3[, tree$tip.label], tree, model.OU.BM)
## [1] -10.71
## attr(,"X0")
## [1]  15.99  18.34 -11.95

Tracing the likelihood calculation using the function PCMLikTrace

Variant 1

traceTable1 <- PCMLikTrace(X[, tree$tip.label], tree, model.OU.BM)
traceTable1[, node:=.I]
setkey(traceTable1, i)

Variant 2

traceTable2 <- PCMLikTrace(
  X[, tree$tip.label], tree, model.OU.BM, metaI = metaI.variant2)
traceTable2[, node:=.I]
setkey(traceTable2, i)

Variant 3

traceTable3 <- PCMLikTrace(X3[, tree$tip.label], tree, model.OU.BM)
traceTable3[, node:=.I]
setkey(traceTable3, i)

A step by step description of the log-likelihood calculation

Step 1: Calculating →ω, Φ and V for each tip or internal node}

Calculating →ω, Φ and V for a node in an OU regime

# OU parameters:
H <- model.OU.BM$`1`$H[,,1]
theta <- model.OU.BM$`1`$Theta[,1]
Sigma <- model.OU.BM$`1`$Sigma_x[,,1] %*% t(model.OU.BM$`1`$Sigma_x[,,1])

# Eigenvalues of H: these can be complex numbers
lambda <- eigen(H)$values

# Matrix of eigenvectors of H: again, these can be complex
P <- eigen(H)$vectors
P_1 <- solve(P)

# vectors of active coordinates:
pc <- PCMInfo(X[, tree$tip.label], tree, model.OU.BM)$pc

# length of the branch leading to tip 1 (2.6):
t1 <- PCMTreeDtNodes(tree)[endNodeLab == "1", endTime - startTime]

# active coordinates for tip 1 and its parent:
k1 <- pc[, match("1", PCMTreeGetLabels(tree))]
k9 <- pc[, match("9", PCMTreeGetLabels(tree))]

# k x k matrix formed from the pairs of lambda-values and t1 (see Eq. 19):
LambdaMat <- matrix(0, 3, 3)
for(i in 1:3) 
  for(j in 1:3) 
    LambdaMat[i,j] <- 1/(lambda[i]+lambda[j])*(1-exp(-(lambda[i]+lambda[j])*t1))

# omega, Phi, V for tip 1:
print(omega1 <- (diag(1, 3, 3)[k1, ] - expm::expm(-H*t1)[k1, ]) %*% theta[])
##        [,1]
## [1,] 0.6151
## [2,] 0.3392
print(Phi1 <- expm::expm(-H*t1)[k1, k9])
##          [,1]    [,2]
## [1,]  0.37818 -0.3461
## [2,] -0.06011  0.1261
print(V1 <- (P %*% (LambdaMat * (P_1 %*% Sigma %*% t(P_1))) %*% t(P))[k1, k1])
##             [,1]        [,2]
## [1,]  2.21599+0i -0.07466+0i
## [2,] -0.07466-0i  0.25787+0i

Calculating →ω, Φ and V for a node in a BM regime

# BM parameter:
Sigma <- model.OU.BM$`2`$Sigma_x[,,1] %*% t(model.OU.BM$`2`$Sigma_x[,,1])

# vectors of active coordinates:
pc <- PCMInfo(X[, tree$tip.label], tree, model.OU.BM)$pc

# length of the branch leading to tip 2 (1.6):
t2 <- PCMTreeDtNodes(tree)[endNodeLab == "2", endTime - startTime]

# active coordinates for tip 1 and its parent:
k2 <- pc[, match("2", PCMTreeGetLabels(tree))]
k6 <- pc[, match("6", PCMTreeGetLabels(tree))]

# omega, Phi, V for tip 1:
print(omega2 <- as.matrix(rep(0, 3)[k2]))
##      [,1]
## [1,]    0
print(Phi2 <- as.matrix(diag(1, 3, 3)[k2, k6]))
##      [,1]
## [1,]    1
## [2,]    0
print(V2 <- as.matrix((t2*Sigma)[k2, k2]))
##      [,1]
## [1,] 2.88

Variant 1

options(digits = 2)
cat(FormatTableAsLatex(
  traceTable1[list(pOrder), list(j, i, t_i, k_i, omega_i, Phi_i, V_i, V_1_i)], 
  type = tableOutputType))









j i ti ki ωi Φi Vi V−1i
6 2 [1.60] [1] [0.00..] [1.00.0.00......] [2.88........] [0.35........]
6 3 [0.80] [13] [0.00.0.00] [1.00.0.00...0.00.1.00] [1.44.0.10...0.10.0.07] [0.76.−1.02...−1.02.15.24]
8 6 [0.70] [13] [0.00.0.00] [1.00.0.00...0.00.1.00] [1.26.0.08...0.08.0.06] [0.87.−1.16...−1.16.17.42]
9 1 [2.60] [13] [0.62.0.34] [0.38.−0.35...−0.06.0.13] [2.22.−0.07...−0.07.0.26] [0.46.0.13...0.13.3.92]
9 8 [0.60] [13] [−0.04.0.49] [0.89.−0.32...−0.17.0.60] [1.00.0.06...0.06.0.21] [1.01.−0.26...−0.26.4.77]
7 4 [1.80] [23] [.−0.860.44] [...0.22−0.21−0.16−0.110.290.24] [....0.37−0.21.−0.210.25] [....5.164.30.4.307.58]
7 5 [0.80] [23] [.−0.780.53] [...0.250.03−0.12−0.170.420.51] [....0.27−0.18.−0.180.22] [....7.786.15.6.159.30]
0 9 [0.90] [13] [0.02.0.53] [0.81−0.61−0.39...−0.170.420.47] [1.35.0.03...0.03.0.23] [0.74.−0.11...−0.11.4.38]
0 7 [1.50] [123] [0.22−0.880.49] [0.64−0.67−0.430.24−0.18−0.16−0.130.340.30] [1.820.45−0.020.450.34−0.20−0.02−0.200.25] [1.29−3.11−2.44−3.1113.0610.44−2.4410.4412.43]
ND 0 ND [123] ND ND ND ND

Variant 2

cat(FormatTableAsLatex(
  traceTable2[list(pOrder), list(j, i, t_i, k_i, omega_i, Phi_i, V_i, V_1_i)], 
  type = tableOutputType))









j i ti ki ωi Φi Vi V−1i
6 2 [1.60] [1] [0.00..] [1.00.0.00......] [2.88........] [0.35........]
6 3 [0.80] [13] [0.00.0.00] [1.00.0.00...0.00.1.00] [1.44.0.10...0.10.0.07] [0.76.−1.02...−1.02.15.24]
8 6 [0.70] [13] [0.00.0.00] [1.00.0.00...0.00.1.00] [1.26.0.08...0.08.0.06] [0.87.−1.16...−1.16.17.42]
9 1 [2.60] [13] [0.62.0.34] [0.38.−0.35...−0.06.0.13] [2.22.−0.07...−0.07.0.26] [0.46.0.13...0.13.3.92]
9 8 [0.60] [13] [−0.04.0.49] [0.89.−0.32...−0.17.0.60] [1.00.0.06...0.06.0.21] [1.01.−0.26...−0.26.4.77]
7 4 [1.80] [23] [.−0.860.44] [...0.22−0.21−0.16−0.110.290.24] [....0.37−0.21.−0.210.25] [....5.164.30.4.307.58]
7 5 [0.80] [23] [.−0.780.53] [...0.250.03−0.12−0.170.420.51] [....0.27−0.18.−0.180.22] [....7.786.15.6.159.30]
0 9 [0.90] [13] [0.02.0.53] [0.81.−0.39...−0.17.0.47] [1.35.0.03...0.03.0.23] [0.74.−0.11...−0.11.4.38]
0 7 [1.50] [123] [0.22−0.880.49] [0.64.−0.430.24.−0.16−0.13.0.30] [1.820.45−0.020.450.34−0.20−0.02−0.200.25] [1.29−3.11−2.44−3.1113.0610.44−2.4410.4412.43]
ND 0 ND [13] ND ND ND ND

Variant 3

options(digits = 2)
cat(FormatTableAsLatex(
  traceTable3[list(pOrder), list(j, i, t_i, k_i, omega_i, Phi_i, V_i, V_1_i)], 
  type = tableOutputType))









j i ti ki ωi Φi Vi V−1i
6 2 [1.60] [1] [0.00..] [1.000.000.00......] [2.88........] [0.35........]
6 3 [0.80] [13] [0.00.0.00] [1.000.000.00...0.000.001.00] [1.44.0.10...0.10.0.07] [0.76.−1.02...−1.02.15.24]
8 6 [0.70] [123] [0.000.000.00] [1.000.000.000.001.000.000.000.001.00] [1.260.350.080.350.240.100.080.100.06] [2.23−7.449.42−7.4440.67−57.879.42−57.8799.76]
9 1 [2.60] [13] [0.62.0.34] [0.38−0.52−0.35...−0.060.170.13] [2.22.−0.07...−0.07.0.26] [0.46.0.13...0.13.3.92]
9 8 [0.60] [123] [−0.04−0.690.49] [0.89−0.51−0.320.230.17−0.10−0.170.400.60] [1.000.190.060.190.24−0.170.06−0.170.21] [2.09−4.60−4.16−4.6019.5016.55−4.1616.5518.82]
7 4 [1.80] [23] [.−0.860.44] [...0.22−0.21−0.16−0.110.290.24] [....0.37−0.21.−0.210.25] [....5.164.30.4.307.58]
7 5 [0.80] [23] [.−0.780.53] [...0.250.03−0.12−0.170.420.51] [....0.27−0.18.−0.180.22] [....7.786.15.6.159.30]
0 9 [0.90] [123] [0.02−0.810.53] [0.81−0.61−0.390.25−0.02−0.13−0.170.420.47] [1.350.290.030.290.28−0.180.03−0.180.23] [1.60−3.66−3.14−3.6615.6212.92−3.1412.9215.08]
0 7 [1.50] [123] [0.22−0.880.49] [0.64−0.67−0.430.24−0.18−0.16−0.130.340.30] [1.820.45−0.020.450.34−0.20−0.02−0.200.25] [1.29−3.11−2.44−3.1113.0610.44−2.4410.4412.43]
ND 0 ND [123] ND ND ND ND

Step 2: Calculating A, →b, C, →d, E and f for tips and internal nodes

# For tip 1. We directly apply Eq. 2, Thm 1:
# We can safely use the real part of V1 (all imaginary parts are 0):
print(V1)
##           [,1]      [,2]
## [1,]  2.216+0i -0.075+0i
## [2,] -0.075-0i  0.258+0i
V1 <- Re(V1)
V1_1 <- solve(V1)

print(A1 <- -0.5*V1_1)
##        [,1]   [,2]
## [1,] -0.228 -0.066
## [2,] -0.066 -1.958
print(E1 <- t(Phi1) %*% V1_1)
##       [,1]  [,2]
## [1,]  0.16 -0.19
## [2,] -0.14  0.45
print(b1 <- V1_1 %*% omega1)
##      [,1]
## [1,] 0.33
## [2,] 1.41
print(C1 <- -0.5 * E1 %*% Phi1)
##        [,1]   [,2]
## [1,] -0.037  0.040
## [2,]  0.040 -0.053
print(d1 <- -E1 %*% omega1)
##        [,1]
## [1,] -0.038
## [2,] -0.065
print(f1 <- -0.5 * (t(omega1) %*% V1_1 %*% omega1 + sum(k1)*log(2*pi) + log(det(V1))))
##      [,1]
## [1,] -1.9

Variant 1

cat(FormatTableAsLatex(
  traceTable1[list(pOrder), list(j, i, k_i, A_i, b_i, C_i, d_i, E_i, f_i)], 
  type = tableOutputType))









j i ki Ai bi Ci di Ei fi
6 2 [1] [−0.17........] [0.00..] [−0.17.−0.00...−0.00.−0.00] [−0.00.−0.00] [0.35.....0.00..] [−1.45]
6 3 [13] [−0.38.0.51...0.51.−7.62] [0.00.0.00] [−0.38.0.51...0.51.−7.62] [0.00.0.00] [0.76.−1.02...−1.02.15.24] [−0.66]
8 6 [13] [−0.44.0.58...0.58.−8.71] [0.00.0.00] [−0.44.0.58...0.58.−8.71] [0.00.0.00] [0.87.−1.16...−1.16.17.42] [−0.52]
9 1 [13] [−0.23.−0.07...−0.07.−1.96] [0.33.1.41] [−0.04.0.04...0.04.−0.05] [−0.04.−0.07] [0.16.−0.19...−0.14.0.45] [−1.89]
9 8 [13] [−0.51.0.13...0.13.−2.39] [−0.17.2.37] [−0.50.0.46...0.46.−0.96] [0.54.−1.48] [0.94.−1.02...−0.48.2.95] [−1.65]
7 4 [23] [....−2.58−2.15.−2.15−3.79] [.−2.54−0.35] [−0.070.050.040.05−0.17−0.140.04−0.14−0.11] [0.53−0.43−0.32] [.0.670.11.0.171.31.0.191.10] [−1.34]
7 5 [23] [....−3.89−3.07.−3.07−4.65] [.−2.850.10] [−0.12−0.000.07−0.00−0.89−0.870.07−0.87−0.89] [0.730.05−0.40] [.0.88−0.06.2.814.07.2.194.00] [−1.21]
0 9 [13] [−0.37.0.06...0.06.−2.19] [−0.04.2.34] [−0.320.360.320.36−0.55−0.550.32−0.55−0.57] [0.43−1.00−1.12] [0.62.−0.83−0.50.1.89−0.34.2.11] [−1.87]
0 7 [123] [−0.641.551.221.55−6.53−5.221.22−5.22−6.22] [1.82−7.04−3.63] [−0.150.230.170.23−0.76−0.580.17−0.58−0.44] [0.061.160.75] [0.40−0.22−0.70−1.133.253.98−0.792.383.09] [−3.47]
ND 0 [123] ND ND ND ND ND ND

Variant 2

cat(FormatTableAsLatex(
  traceTable2[list(pOrder), list(j, i, k_i, A_i, b_i, C_i, d_i, E_i, f_i)], 
  type = tableOutputType))









j i ki Ai bi Ci di Ei fi
6 2 [1] [−0.17........] [0.00..] [−0.17.−0.00...−0.00.−0.00] [−0.00.−0.00] [0.35.....0.00..] [−1.45]
6 3 [13] [−0.38.0.51...0.51.−7.62] [0.00.0.00] [−0.38.0.51...0.51.−7.62] [0.00.0.00] [0.76.−1.02...−1.02.15.24] [−0.66]
8 6 [13] [−0.44.0.58...0.58.−8.71] [0.00.0.00] [−0.44.0.58...0.58.−8.71] [0.00.0.00] [0.87.−1.16...−1.16.17.42] [−0.52]
9 1 [13] [−0.23.−0.07...−0.07.−1.96] [0.33.1.41] [−0.04.0.04...0.04.−0.05] [−0.04.−0.07] [0.16.−0.19...−0.14.0.45] [−1.89]
9 8 [13] [−0.51.0.13...0.13.−2.39] [−0.17.2.37] [−0.50.0.46...0.46.−0.96] [0.54.−1.48] [0.94.−1.02...−0.48.2.95] [−1.65]
7 4 [23] [....−2.58−2.15.−2.15−3.79] [.−2.54−0.35] [−0.070.050.040.05−0.17−0.140.04−0.14−0.11] [0.53−0.43−0.32] [.0.670.11.0.171.31.0.191.10] [−1.34]
7 5 [23] [....−3.89−3.07.−3.07−4.65] [.−2.850.10] [−0.12−0.000.07−0.00−0.89−0.870.07−0.87−0.89] [0.730.05−0.40] [.0.88−0.06.2.814.07.2.194.00] [−1.21]
0 9 [13] [−0.37.0.06...0.06.−2.19] [−0.04.2.34] [−0.32.0.32...0.32.−0.57] [0.43.−1.12] [0.62.−0.83...−0.34.2.11] [−1.87]
0 7 [123] [−0.641.551.221.55−6.53−5.221.22−5.22−6.22] [1.82−7.04−3.63] [−0.15.0.17...0.17.−0.44] [0.06.0.75] [0.40−0.22−0.70...−0.792.383.09] [−3.47]
ND 0 [13] ND ND ND ND ND ND

Variant 3

cat(FormatTableAsLatex(
  traceTable3[list(pOrder), list(j, i, k_i, A_i, b_i, C_i, d_i, E_i, f_i)], 
  type = tableOutputType))









j i ki Ai bi Ci di Ei fi
6 2 [1] [−0.17........] [0.00..] [−0.17−0.00−0.00−0.00−0.00−0.00−0.00−0.00−0.00] [−0.00−0.00−0.00] [0.35..0.00..0.00..] [−1.45]
6 3 [13] [−0.38.0.51...0.51.−7.62] [0.00.0.00] [−0.38−0.000.51−0.00−0.00−0.000.51−0.00−7.62] [0.000.000.00] [0.76.−1.020.00.0.00−1.02.15.24] [−0.66]
8 6 [123] [−1.123.72−4.713.72−20.3428.94−4.7128.94−49.88] [0.000.000.00] [−1.123.72−4.713.72−20.3428.94−4.7128.94−49.88] [0.000.000.00] [2.23−7.449.42−7.4440.67−57.879.42−57.8799.76] [0.41]
9 1 [13] [−0.23.−0.07...−0.07.−1.96] [0.33.1.41] [−0.040.060.040.06−0.11−0.080.04−0.08−0.05] [−0.04−0.07−0.07] [0.16.−0.19−0.22.0.61−0.14.0.45] [−1.89]
9 8 [123] [−1.052.302.082.30−9.75−8.282.08−8.28−9.41] [1.04−5.12−1.98] [−0.651.191.041.19−4.36−3.661.04−3.66−3.26] [−0.082.171.01] [1.50−2.38−3.04−3.4912.1712.36−2.719.4610.98] [−1.75]
7 4 [23] [....−2.58−2.15.−2.15−3.79] [.−2.54−0.35] [−0.070.050.040.05−0.17−0.140.04−0.14−0.11] [0.53−0.43−0.32] [.0.670.11.0.171.31.0.191.10] [−1.34]
7 5 [23] [....−3.89−3.07.−3.07−4.65] [.−2.850.10] [−0.12−0.000.07−0.00−0.89−0.870.07−0.87−0.89] [0.730.05−0.40] [.0.88−0.06.2.814.07.2.194.00] [−1.21]
0 9 [123] [−0.801.831.571.83−7.81−6.461.57−6.46−7.54] [1.34−5.90−2.54] [−0.360.640.520.64−2.26−1.830.52−1.83−1.53] [−0.011.760.96] [0.89−1.17−1.80−2.227.317.94−1.635.506.66] [−2.53]
0 7 [123] [−0.641.551.221.55−6.53−5.221.22−5.22−6.22] [1.82−7.04−3.63] [−0.150.230.170.23−0.76−0.580.17−0.58−0.44] [0.061.160.75] [0.40−0.22−0.70−1.133.253.98−0.792.383.09] [−3.47]
ND 0 [123] ND ND ND ND ND ND

Step 3: Calculating L, →m, r for the internal nodes and the root

# For tip 2 with parent node 6, we use the following terms stored in Table S5:
A2 <- matrix(-0.17)
b2 <- 0.0
C2 <- rbind(c(-0.17, 0), 
            c(0, 0))
d2 <- c(0.0, 0.0)
E2 <- matrix(c(0.35, 0), nrow = 2, ncol = 1)
f2 <- -1.45
k2 <- 1

# Now we apply Eq. S3:
print(L62 <- C2)
##       [,1] [,2]
## [1,] -0.17    0
## [2,]  0.00    0
print(m62 <- d2 + E2 %*% X[k2, "2", drop = FALSE])
##          2
## [1,] 0.035
## [2,] 0.000
print(r62 <- t(X[k2, "2", drop = FALSE]) %*% A2 %*% X[k2, "2", drop = FALSE] + 
        t(X[k2, "2", drop = FALSE]) %*% b2 + f2)
##      2
## 2 -1.5
# For tip 3 with parent node 6, applying Eq. S3, we obtain (see Table S8):
L63 <- rbind(c(-0.38, 0.51),
             c(0.51, -7.62))
m63 <- c(-1.07, 18.09)
r63 <- -11.41

# Now, we sum the terms L6i, m6i and r6i over all daughters of 6 (i) to obtain:
print(L6 <- L62 + L63)
##       [,1]  [,2]
## [1,] -0.55  0.51
## [2,]  0.51 -7.62
print(m6 <- m62 + m63)
##       2
## [1,] -1
## [2,] 18
print(r6 <- r62 + r63)
##     2
## 2 -13
# Using Eq. S3, we obtain L86, m86, r86, using the values for A,b,C,d,E,f in Table S5:
A6 <- rbind(c(-0.44, 0.58),
            c(0.58, -8.71))
b6 <- c(0.0, 0.0)
C6 <- rbind(c(-0.44, 0.58),
            c(0.58, -8.71))
d6 <- c(0.0, 0.0)
E6 <- rbind(c(0.87, -1.16),
            c(-1.16, 17.42))
f6 <- -0.52
k6 <- c(1, 3)

print(L86 <- C6 - (1/4)*E6 %*% solve(A6 + L6) %*% t(E6))
##       [,1]  [,2]
## [1,] -0.25  0.27
## [2,]  0.27 -4.06
print(m86 <- d6 - (1/2)*E6 %*% solve(A6 + L6) %*% (b6+m6))
##          2
## [1,] -0.57
## [2,]  9.65
print(r86 <- f6+r6+(length(k6)/2)*log(2*pi) - 
        (1/2)*log(det(-2*(A6+L6))) -
        (1/4)*t(b6+m6) %*% solve(A6+L6) %*% (b6+m6))
##      2
## 2 -8.6
# Because 8 is a singleton node, we immediately obtain L8, m8, r8:
L8 <- L86; m8 <- m86; r8 <- r86;

Variant 1

options(digits = 3)
cat(FormatTableAsLatex(
  traceTable1[
    list(pOrder), 
    list(
      j, i, X_i, k_i, 
      L_i, m_i, r_i,
      `L_{ji}`, `m_{ji}`, `r_{ji}`)], 
  
  type = tableOutputType))









j i Xi ki Li mi ri Lji mji rji
6 2 [0.1NaNNA] [1] ND ND ND [−0.174.−0.000...−0.000.−0.000] [0.035.0.000] [−1.450]
6 3 [0.2NaN1.2] [13] ND ND ND [−0.381.0.508...0.508.−7.622] [−1.067.18.089] [−11.405]
8 6 [NANaNNA] [13] [−0.555.0.508...0.508.−7.622] [−1.032.18.089] [−12.855] [−0.243.0.271...0.271.−4.065] [−0.568.9.648] [−8.571]
9 1 [0.3NaN1.4] [13] ND ND ND [−0.037.0.040...0.040.−0.053] [−0.249.0.520] [−3.735]
9 8 [NANaNNA] [13] [−0.243.0.271...0.271.−4.065] [−0.568.9.648] [−8.571] [−0.195.0.250...0.250.−0.594] [−0.402.1.267] [−4.248]
7 4 [NA0.20.2] [23] ND ND ND [−0.0680.0540.0400.054−0.174−0.1410.040−0.141−0.114] [0.683−0.138−0.066] [−2.349]
7 5 [NA1.20.4] [23] ND ND ND [−0.115−0.0010.070−0.001−0.893−0.8690.070−0.869−0.890] [1.7625.0433.834] [−13.880]
0 9 [NANaNNA] [13] [−0.232.0.291...0.291.−0.647] [−0.651.1.786] [−7.983] [−0.1400.1620.1430.162−0.200−0.1830.143−0.183−0.170] [−0.2620.4220.429] [−7.430]
0 7 [NANANA] [123] [−0.1830.0530.1100.053−1.066−1.0100.110−1.010−1.004] [2.4444.9063.769] [−16.230] [−0.0520.0520.0350.052−0.113−0.0820.035−0.082−0.060] [1.222−0.396−0.174] [−10.947]
ND 0 [NANANA] [123] [−0.1920.2140.1780.214−0.313−0.2650.178−0.265−0.230] [0.9600.0260.255] [−18.377] ND ND ND

Variant 2

cat(FormatTableAsLatex(
  traceTable2[
    list(pOrder), 
    list(
      j, i, X_i, k_i, 
      L_i, m_i, r_i,
      `L_{ji}`, `m_{ji}`, `r_{ji}`)], 
  
  type = tableOutputType))









j i Xi ki Li mi ri Lji mji rji
6 2 [0.1NaNNA] [1] ND ND ND [−0.174.−0.000...−0.000.−0.000] [0.035.0.000] [−1.450]
6 3 [0.2NaN1.2] [13] ND ND ND [−0.381.0.508...0.508.−7.622] [−1.067.18.089] [−11.405]
8 6 [NANaNNA] [13] \begin{bmatrix}{} -0.555 & . & 0.508 \\ . & . & . \\ 0.508 & . & -7.622 \\ \end{bmatrix} \begin{bmatrix}{} -1.032 \\ . \\ 18.089 \\ \end{bmatrix} \begin{bmatrix}{} -12.855 \\ \end{bmatrix} \begin{bmatrix}{} -0.243 & . & 0.271 \\ . & . & . \\ 0.271 & . & -4.065 \\ \end{bmatrix} \begin{bmatrix}{} -0.568 \\ . \\ 9.648 \\ \end{bmatrix} \begin{bmatrix}{} -8.571 \\ \end{bmatrix}
9 1 \begin{bmatrix}{} 0.3 \\ NaN \\ 1.4 \\ \end{bmatrix} \begin{bmatrix}{} 1 \\ 3 \\ \end{bmatrix} ND ND ND \begin{bmatrix}{} -0.037 & . & 0.040 \\ . & . & . \\ 0.040 & . & -0.053 \\ \end{bmatrix} \begin{bmatrix}{} -0.249 \\ . \\ 0.520 \\ \end{bmatrix} \begin{bmatrix}{} -3.735 \\ \end{bmatrix}
9 8 \begin{bmatrix}{} NA \\ NaN \\ NA \\ \end{bmatrix} \begin{bmatrix}{} 1 \\ 3 \\ \end{bmatrix} \begin{bmatrix}{} -0.243 & . & 0.271 \\ . & . & . \\ 0.271 & . & -4.065 \\ \end{bmatrix} \begin{bmatrix}{} -0.568 \\ . \\ 9.648 \\ \end{bmatrix} \begin{bmatrix}{} -8.571 \\ \end{bmatrix} \begin{bmatrix}{} -0.195 & . & 0.250 \\ . & . & . \\ 0.250 & . & -0.594 \\ \end{bmatrix} \begin{bmatrix}{} -0.402 \\ . \\ 1.267 \\ \end{bmatrix} \begin{bmatrix}{} -4.248 \\ \end{bmatrix}
7 4 \begin{bmatrix}{} NA \\ 0.2 \\ 0.2 \\ \end{bmatrix} \begin{bmatrix}{} 2 \\ 3 \\ \end{bmatrix} ND ND ND \begin{bmatrix}{} -0.068 & 0.054 & 0.040 \\ 0.054 & -0.174 & -0.141 \\ 0.040 & -0.141 & -0.114 \\ \end{bmatrix} \begin{bmatrix}{} 0.683 \\ -0.138 \\ -0.066 \\ \end{bmatrix} \begin{bmatrix}{} -2.349 \\ \end{bmatrix}
7 5 \begin{bmatrix}{} NA \\ 1.2 \\ 0.4 \\ \end{bmatrix} \begin{bmatrix}{} 2 \\ 3 \\ \end{bmatrix} ND ND ND \begin{bmatrix}{} -0.115 & -0.001 & 0.070 \\ -0.001 & -0.893 & -0.869 \\ 0.070 & -0.869 & -0.890 \\ \end{bmatrix} \begin{bmatrix}{} 1.762 \\ 5.043 \\ 3.834 \\ \end{bmatrix} \begin{bmatrix}{} -13.880 \\ \end{bmatrix}
0 9 \begin{bmatrix}{} NA \\ NaN \\ NA \\ \end{bmatrix} \begin{bmatrix}{} 1 \\ 3 \\ \end{bmatrix} \begin{bmatrix}{} -0.232 & . & 0.291 \\ . & . & . \\ 0.291 & . & -0.647 \\ \end{bmatrix} \begin{bmatrix}{} -0.651 \\ . \\ 1.786 \\ \end{bmatrix} \begin{bmatrix}{} -7.983 \\ \end{bmatrix} \begin{bmatrix}{} -0.140 & . & 0.143 \\ . & . & . \\ 0.143 & . & -0.170 \\ \end{bmatrix} \begin{bmatrix}{} -0.262 \\ . \\ 0.429 \\ \end{bmatrix} \begin{bmatrix}{} -7.430 \\ \end{bmatrix}
0 7 \begin{bmatrix}{} NA \\ NA \\ NA \\ \end{bmatrix} \begin{bmatrix}{} 1 \\ 2 \\ 3 \\ \end{bmatrix} \begin{bmatrix}{} -0.183 & 0.053 & 0.110 \\ 0.053 & -1.066 & -1.010 \\ 0.110 & -1.010 & -1.004 \\ \end{bmatrix} \begin{bmatrix}{} 2.444 \\ 4.906 \\ 3.769 \\ \end{bmatrix} \begin{bmatrix}{} -16.230 \\ \end{bmatrix} \begin{bmatrix}{} -0.052 & 0.053 & 0.035 \\ 0.053 & -1.066 & -1.010 \\ 0.035 & -1.010 & -0.060 \\ \end{bmatrix} \begin{bmatrix}{} 1.222 \\ 4.906 \\ -0.174 \\ \end{bmatrix} \begin{bmatrix}{} -10.947 \\ \end{bmatrix}
ND 0 \begin{bmatrix}{} NA \\ NaN \\ NA \\ \end{bmatrix} \begin{bmatrix}{} 1 \\ 3 \\ \end{bmatrix} \begin{bmatrix}{} -0.192 & . & 0.178 \\ . & . & . \\ 0.178 & . & -0.230 \\ \end{bmatrix} \begin{bmatrix}{} 0.960 \\ . \\ 0.255 \\ \end{bmatrix} \begin{bmatrix}{} -18.377 \\ \end{bmatrix} ND ND ND

Variant 3

cat(FormatTableAsLatex(
  traceTable3[
    list(pOrder), 
    list(
      j, i, X_i, k_i, 
      L_i, m_i, r_i, 
      `L_{ji}`, `m_{ji}`, `r_{ji}`)], 
  
  type = tableOutputType))









j i X_i k_i L_i m_i r_i L_{ji} m_{ji} r_{ji}
6 2 \begin{bmatrix}{} 0.1 \\ NA \\ NA \\ \end{bmatrix} \begin{bmatrix}{} 1 \\ \end{bmatrix} ND ND ND \begin{bmatrix}{} -0.174 & -0.000 & -0.000 \\ -0.000 & -0.000 & -0.000 \\ -0.000 & -0.000 & -0.000 \\ \end{bmatrix} \begin{bmatrix}{} 0.035 \\ 0.000 \\ 0.000 \\ \end{bmatrix} \begin{bmatrix}{} -1.450 \\ \end{bmatrix}
6 3 \begin{bmatrix}{} 0.2 \\ NA \\ 1.2 \\ \end{bmatrix} \begin{bmatrix}{} 1 \\ 3 \\ \end{bmatrix} ND ND ND \begin{bmatrix}{} -0.381 & -0.000 & 0.508 \\ -0.000 & -0.000 & -0.000 \\ 0.508 & -0.000 & -7.622 \\ \end{bmatrix} \begin{bmatrix}{} -1.067 \\ 0.000 \\ 18.089 \\ \end{bmatrix} \begin{bmatrix}{} -11.405 \\ \end{bmatrix}
8 6 \begin{bmatrix}{} NA \\ NA \\ NA \\ \end{bmatrix} \begin{bmatrix}{} 1 \\ 2 \\ 3 \\ \end{bmatrix} \begin{bmatrix}{} -0.555 & 0.000 & 0.508 \\ 0.000 & 0.000 & 0.000 \\ 0.508 & 0.000 & -7.622 \\ \end{bmatrix} \begin{bmatrix}{} -1.032 \\ 0.000 \\ 18.089 \\ \end{bmatrix} \begin{bmatrix}{} -12.855 \\ \end{bmatrix} \begin{bmatrix}{} -0.243 & -0.000 & 0.271 \\ -0.000 & 0.000 & -0.000 \\ 0.271 & -0.000 & -4.065 \\ \end{bmatrix} \begin{bmatrix}{} -0.568 \\ -0.000 \\ 9.648 \\ \end{bmatrix} \begin{bmatrix}{} -8.571 \\ \end{bmatrix}
9 1 \begin{bmatrix}{} 0.3 \\ NA \\ 1.4 \\ \end{bmatrix} \begin{bmatrix}{} 1 \\ 3 \\ \end{bmatrix} ND ND ND \begin{bmatrix}{} -0.037 & 0.059 & 0.040 \\ 0.059 & -0.109 & -0.076 \\ 0.040 & -0.076 & -0.053 \\ \end{bmatrix} \begin{bmatrix}{} -0.249 \\ 0.711 \\ 0.520 \\ \end{bmatrix} \begin{bmatrix}{} -3.735 \\ \end{bmatrix}
9 8 \begin{bmatrix}{} NA \\ NA \\ NA \\ \end{bmatrix} \begin{bmatrix}{} 1 \\ 2 \\ 3 \\ \end{bmatrix} \begin{bmatrix}{} -0.243 & -0.000 & 0.271 \\ -0.000 & 0.000 & -0.000 \\ 0.271 & -0.000 & -4.065 \\ \end{bmatrix} \begin{bmatrix}{} -0.568 \\ -0.000 \\ 9.648 \\ \end{bmatrix} \begin{bmatrix}{} -8.571 \\ \end{bmatrix} \begin{bmatrix}{} -0.195 & 0.213 & 0.250 \\ 0.213 & -0.318 & -0.426 \\ 0.250 & -0.426 & -0.594 \\ \end{bmatrix} \begin{bmatrix}{} -0.402 \\ 0.860 \\ 1.267 \\ \end{bmatrix} \begin{bmatrix}{} -4.248 \\ \end{bmatrix}
7 4 \begin{bmatrix}{} NA \\ 0.2 \\ 0.2 \\ \end{bmatrix} \begin{bmatrix}{} 2 \\ 3 \\ \end{bmatrix} ND ND ND \begin{bmatrix}{} -0.068 & 0.054 & 0.040 \\ 0.054 & -0.174 & -0.141 \\ 0.040 & -0.141 & -0.114 \\ \end{bmatrix} \begin{bmatrix}{} 0.683 \\ -0.138 \\ -0.066 \\ \end{bmatrix} \begin{bmatrix}{} -2.349 \\ \end{bmatrix}
7 5 \begin{bmatrix}{} NA \\ 1.2 \\ 0.4 \\ \end{bmatrix} \begin{bmatrix}{} 2 \\ 3 \\ \end{bmatrix} ND ND ND \begin{bmatrix}{} -0.115 & -0.001 & 0.070 \\ -0.001 & -0.893 & -0.869 \\ 0.070 & -0.869 & -0.890 \\ \end{bmatrix} \begin{bmatrix}{} 1.762 \\ 5.043 \\ 3.834 \\ \end{bmatrix} \begin{bmatrix}{} -13.880 \\ \end{bmatrix}
0 9 \begin{bmatrix}{} NA \\ NA \\ NA \\ \end{bmatrix} \begin{bmatrix}{} 1 \\ 2 \\ 3 \\ \end{bmatrix} \begin{bmatrix}{} -0.232 & 0.272 & 0.291 \\ 0.272 & -0.427 & -0.502 \\ 0.291 & -0.502 & -0.647 \\ \end{bmatrix} \begin{bmatrix}{} -0.651 \\ 1.571 \\ 1.786 \\ \end{bmatrix} \begin{bmatrix}{} -7.983 \\ \end{bmatrix} \begin{bmatrix}{} -0.089 & 0.140 & 0.106 \\ 0.140 & -0.258 & -0.204 \\ 0.106 & -0.204 & -0.163 \\ \end{bmatrix} \begin{bmatrix}{} -0.394 \\ 1.044 \\ 0.833 \\ \end{bmatrix} \begin{bmatrix}{} -8.380 \\ \end{bmatrix}
0 7 \begin{bmatrix}{} NA \\ NA \\ NA \\ \end{bmatrix} \begin{bmatrix}{} 1 \\ 2 \\ 3 \\ \end{bmatrix} \begin{bmatrix}{} -0.183 & 0.053 & 0.110 \\ 0.053 & -1.066 & -1.010 \\ 0.110 & -1.010 & -1.004 \\ \end{bmatrix} \begin{bmatrix}{} 2.444 \\ 4.906 \\ 3.769 \\ \end{bmatrix} \begin{bmatrix}{} -16.230 \\ \end{bmatrix} \begin{bmatrix}{} -0.052 & 0.052 & 0.035 \\ 0.052 & -0.113 & -0.082 \\ 0.035 & -0.082 & -0.060 \\ \end{bmatrix} \begin{bmatrix}{} 1.222 \\ -0.396 \\ -0.174 \\ \end{bmatrix} \begin{bmatrix}{} -10.947 \\ \end{bmatrix}
ND 0 \begin{bmatrix}{} NA \\ NA \\ NA \\ \end{bmatrix} \begin{bmatrix}{} 1 \\ 2 \\ 3 \\ \end{bmatrix} \begin{bmatrix}{} -0.141 & 0.192 & 0.141 \\ 0.192 & -0.372 & -0.286 \\ 0.141 & -0.286 & -0.223 \\ \end{bmatrix} \begin{bmatrix}{} 0.827 \\ 0.648 \\ 0.659 \\ \end{bmatrix} \begin{bmatrix}{} -19.328 \\ \end{bmatrix} ND ND ND

Step 4: Calculating the log-likelihood value using \mathbf{L}_{0}, \vec{m}_{0}, and r_{0}.

# Variant 1.
# Copy the values of L0, m0 and r0 from Table S8:
L0 <- rbind(c(-0.192, 0.214, 0.178),
            c(0.214, -0.313, -0.265),
            c(0.178, -0.265, -0.230))
m0 <- c(0.96, 0.026, 0.255)
r0 <- -18.377

# Use Eq. S2 to estimate the optimal X0:
print(t(x0Hat <- -0.5*solve(L0) %*% m0))
##      [,1]  [,2] [,3]
## [1,] 9.64 -6.25 15.2
# Use Eq. S1 to calculate the log-likelihood:
print(ll0 <- t(x0Hat) %*% L0 %*% x0Hat + t(x0Hat) %*% m0 + r0)
##       [,1]
## [1,] -11.9
# Variant 2.
# Copy the values of L0, m0 and r0 from Table S9:
L0 <- rbind(c(-0.192, 0.178),
            c(0.178, -0.230))
m0 <- c(0.96, 0.255)
r0 <- -18.377

# Use Eq. S2 to estimate the optimal X0:
print(t(x0Hat <- -0.5*solve(L0) %*% m0))
##      [,1] [,2]
## [1,] 10.7 8.81
# Use Eq. S1 to calculate the log-likelihood:
print(ll0 <- t(x0Hat) %*% L0 %*% x0Hat + t(x0Hat) %*% m0 + r0)
##       [,1]
## [1,] -12.1
# Variant 3.
# The function PCMLikTrace generates a data.table with the values of 
# omega, Phi, V, A, b, C, d, E, f, L, m, r. 
traceTable3 <- PCMLikTrace(X3[, tree$tip.label], tree, model.OU.BM)
# The column i corresponds to the node label in the tree as depicted on Fig. 1:
setkey(traceTable3, i)

options(digits = 4)
# Variant 3.
# Copy the values of L0, m0 and r0 from the traceTable object (these values have
# the maximal double floating point precision):
print(L0 <- traceTable3[list("0")][["L_i"]][[1]])
##         [,1]    [,2]    [,3]
## [1,] -0.1408  0.1919  0.1407
## [2,]  0.1919 -0.3715 -0.2862
## [3,]  0.1407 -0.2862 -0.2233
print(m0 <- traceTable3[list("0")][["m_i"]][[1]])
## [1] 0.8273 0.6484 0.6590
print(r0 <- traceTable3[list("0")][["r_i"]][[1]])
## [1] -19.33
# Notice the exact match with the values for variant 3 reported in Fig. S5:
print(t(x0Hat <- -0.5*solve(L0) %*% m0))
##       [,1]  [,2]   [,3]
## [1,] 15.99 18.34 -11.95
print(ll0 <- t(x0Hat) %*% L0 %*% x0Hat + t(x0Hat) %*% m0 + r0)
##        [,1]
## [1,] -10.71

References

Mitov, Venelin, Krzysztof Bartoszek, Georgios Asimomitis, and Tanja Stadler. 2019. “Fast likelihood calculation for multivariate Gaussian phylogenetic models with shifts.” Theor. Popul. Biol., December. https://doi.org/10.1016/j.tpb.2019.11.005.