---
title: "Tracing the likelihood calculation of a Gaussian model"
author: "Venelin Mitov"
date: '`r Sys.Date()`'
output: 
  rmarkdown::html_vignette:
    fig_caption: yes
    toc: true
    toc_depth: 3
vignette: >
  %\VignetteIndexEntry{Tracing the likelihood calculation of a Gaussian model}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
bibliography: REFERENCES.bib
---

<!--
# Copyright 2016-2021 Venelin Mitov
#
# This file is part of PCMBase.
#
# PCMBase is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# PCMBase is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with PCMBase.  If not, see <http://www.gnu.org/licenses/>.
-->

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(ape)
library(PCMBase)
library(data.table)
library(xtable)

FLAGSuggestsAvailable <- PCMBase::RequireSuggestedPackages()

options(digits = 2)
# specify either 'html' or 'latex'
tableOutputType <- 'html' # 'latex' is used for generating the S-tables in the ms.
```

This tutorial shows how the pruning likelihood calculation algorithm works for $\mathcal{G}_{LInv}$ 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:2018fl].

# Example tree and data
```{r treeAndData}
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)
```

```{r PlotTreeAndData, include=FALSE, eval=FALSE}
library(tikzDevice); library(ggplot2); library(data.table); 
# 4. Plotting the tree, the data and the active coordinate vectors:
tipValueLabels <- data.table(
  node = seq_len(PCMTreeNumTips(tree)),
  valueLabel = paste0(
    "$\\vec{x}_{", tree$tip.label, "}=(", apply(X[, tree$tip.label], 2, toString), ")^T$"),
  parse = TRUE)

# Determine the active coordinates for X:
k_i <- PCMPresentCoordinates(X[, tree$tip.label], tree, NULL)


dtNodes <- PCMTreeDtNodes(tree)
dtNodes[, kLabel:=paste0(
  "$\\vec{k}_{", endNodeLab, "}=(", sapply(endNode, function(i) toString(which(k_i[,i]))), ")^T$")]
dtNodes[, kLabel2:=kLabel]
dtNodes[endNodeLab == "8", kLabel:=NA]
dtNodes[endNodeLab != "8", kLabel2:=NA]
dtNodes[, tLabel:=paste0("$t_{", endNodeLab, "}=", endTime-startTime, "$")]
dtNodes[endNodeLab == "0", tLabel:=NA]

tikz(file = "TreeMGPMExample.tex", width = 8, height = 5)
palette <- PCMColorPalette(2, names = c("1", "2"), colors = c("black", "orange"))

plTree <- PCMTreePlot(tree, palette = palette, size=2)

# Plot the tree with branches colored according to the regimes.
# The following code works correctly only if the ggtree package is installed, 
# which is not on CRAN.
if(requireNamespace("ggtree")) {
  
  plTree <- plTree + 
    ggtree::geom_nodelab(geom = "label", color = "red") + 
    ggtree::geom_tiplab(geom = "label", color = "black")
  
  plTree <- plTree %<+% tipValueLabels %<+% dtNodes[, list(node = endNode, kLabel, kLabel2, tLabel)]
  
  plTree <- plTree +
    ggtree::geom_tiplab(geom = "text", aes(label = valueLabel), color = "black", hjust = -0.2, vjust = -1.1) +
    ggtree::geom_tiplab(geom = "text", aes(label = kLabel), color = "black", hjust = -0.4, vjust = 1.1) +
    ggtree::geom_nodelab(geom = "text", aes(label = kLabel), color = "red", hjust = -0.2, vjust = 0.6) +
    ggtree::geom_nodelab(geom = "text", aes(label = kLabel2), color = "red", hjust = 0.4, vjust = 2.8) +
    geom_text(aes(x = branch, label = tLabel), vjust = -0.8, color = "black") +
    scale_x_continuous(limits = c(0, 5.2)) + scale_y_continuous(limits = c(0.8, 5.2))
} 

plTree

dev.off()
```

![A tree with five tips and two evolutionary regimes](Fig1.pdf){height="500px" width="100%"}

# A mixed Gaussian model
```{r, results='asis'}
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)
```

```{r add-to-PCMBaseTestObjects, include = FALSE, eval=FALSE}
# add these objects to the PCMBaseTestObjects (needed for the coding examples).
PCMBaseTestObjects[["tree"]] <- tree
PCMBaseTestObjects[["X"]] <- X[, tree$tip.label]
PCMBaseTestObjects[["model.OU.BM"]] <- model.OU.BM

usethis::use_data(PCMBaseTestObjects, overwrite = TRUE)
```

# Likelihood values for the three example variants
```{r, echo = TRUE}
options(digits = 4)
# Variant 1: 
PCMLik(X[, tree$tip.label], tree, model.OU.BM)

# 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
# And the log-likelihood value is:
PCMLik(X[, tree$tip.label], tree, model.OU.BM, metaI = metaI.variant2)

# 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)
```

# Tracing the likelihood calculation using the function `PCMLikTrace`

## Variant 1
```{r traceTable1R}
traceTable1 <- PCMLikTrace(X[, tree$tip.label], tree, model.OU.BM)
traceTable1[, node:=.I]
setkey(traceTable1, i)
```

## Variant 2
```{r traceTable2R}
traceTable2 <- PCMLikTrace(
  X[, tree$tip.label], tree, model.OU.BM, metaI = metaI.variant2)
traceTable2[, node:=.I]
setkey(traceTable2, i)
```

## Variant 3
```{r traceTable3R}
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 $\vec{\omega}$, $\mathbf{\Phi}$ and $\mathbf{V}$ for each tip or internal node}

### Calculating $\vec{\omega}$, $\mathbf{\Phi}$ and $\mathbf{V}$ for a node in an OU regime
```{r omegaPhiVOU}
# 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[])
print(Phi1 <- expm::expm(-H*t1)[k1, k9])
print(V1 <- (P %*% (LambdaMat * (P_1 %*% Sigma %*% t(P_1))) %*% t(P))[k1, k1])
```

### Calculating $\vec{\omega}$, $\mathbf{\Phi}$ and $\mathbf{V}$ for a node in a BM regime
```{r omegaPhiVBM}
# 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]))
print(Phi2 <- as.matrix(diag(1, 3, 3)[k2, k6]))
print(V2 <- as.matrix((t2*Sigma)[k2, k2]))
```

### Variant 1
```{r omegaPhiV1R, results='asis'}
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))
```

### Variant 2
```{r omegaPhiV2R, results='asis'}
cat(FormatTableAsLatex(
  traceTable2[list(pOrder), list(j, i, t_i, k_i, omega_i, Phi_i, V_i, V_1_i)], 
  type = tableOutputType))
```

### Variant 3
```{r omegaPhiV3R, results='asis'}
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))
```

## Step 2: Calculating $\mathbf{A}$, $\vec{b}$, $\mathbf{C}$, $\vec{d}$, $\mathbf{E}$ and $f$ for tips and internal nodes

```{r AbCdEf}
# 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)
V1 <- Re(V1)
V1_1 <- solve(V1)

print(A1 <- -0.5*V1_1)
print(E1 <- t(Phi1) %*% V1_1)
print(b1 <- V1_1 %*% omega1)
print(C1 <- -0.5 * E1 %*% Phi1)
print(d1 <- -E1 %*% omega1)
print(f1 <- -0.5 * (t(omega1) %*% V1_1 %*% omega1 + sum(k1)*log(2*pi) + log(det(V1))))
```

### Variant 1
```{r AbCdEf1R, results='asis'}
cat(FormatTableAsLatex(
  traceTable1[list(pOrder), list(j, i, k_i, A_i, b_i, C_i, d_i, E_i, f_i)], 
  type = tableOutputType))
```

### Variant 2
```{r AbCdEf2R, results='asis'}
cat(FormatTableAsLatex(
  traceTable2[list(pOrder), list(j, i, k_i, A_i, b_i, C_i, d_i, E_i, f_i)], 
  type = tableOutputType))
```

### Variant 3
```{r AbCdEf3R, results='asis'}
cat(FormatTableAsLatex(
  traceTable3[list(pOrder), list(j, i, k_i, A_i, b_i, C_i, d_i, E_i, f_i)], 
  type = tableOutputType))
```

## Step 3: Calculating $\mathbf{L}$, $\vec{m}$, $r$ for the internal nodes and the root
```{r LmrVariant1}
# 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)
print(m62 <- d2 + E2 %*% X[k2, "2", drop = FALSE])
print(r62 <- t(X[k2, "2", drop = FALSE]) %*% A2 %*% X[k2, "2", drop = FALSE] + 
        t(X[k2, "2", drop = FALSE]) %*% b2 + f2)

# 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)
print(m6 <- m62 + m63)
print(r6 <- r62 + r63)

# 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))
print(m86 <- d6 - (1/2)*E6 %*% solve(A6 + L6) %*% (b6+m6))
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))

# Because 8 is a singleton node, we immediately obtain L8, m8, r8:
L8 <- L86; m8 <- m86; r8 <- r86;
```

### Variant 1
```{r Lmr1R, results='asis'}
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))
```

### Variant 2
```{r Lmr2R, results='asis'}
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))
```

### Variant 3
```{r Lmr3R, results='asis'}
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))
```

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

```{r}
# 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))
# Use Eq. S1 to calculate the log-likelihood:
print(ll0 <- t(x0Hat) %*% L0 %*% x0Hat + t(x0Hat) %*% m0 + r0)

# 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))

# Use Eq. S1 to calculate the log-likelihood:
print(ll0 <- t(x0Hat) %*% L0 %*% x0Hat + t(x0Hat) %*% m0 + r0)

# 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]])
print(m0 <- traceTable3[list("0")][["m_i"]][[1]])
print(r0 <- traceTable3[list("0")][["r_i"]][[1]])

# Notice the exact match with the values for variant 3 reported in Fig. S5:
print(t(x0Hat <- -0.5*solve(L0) %*% m0))
print(ll0 <- t(x0Hat) %*% L0 %*% x0Hat + t(x0Hat) %*% m0 + r0)
```

# References