treeppl_x_revbayes.qmd

---
title: "Comparison of inference under the host repertoire model in TreePPL vs RevBayes"
author: "Mariana Braga"
date: last-modified
date-format: long
format: 
  html:
    code-fold: true
    code-summary: "Show the code"
editor: source
editor_options: 
  chunk_output_type: console
execute:
  warning: false
---

Packages

```{r}
#| label: "packages"
#| message: FALSE

#install_github("maribraga/treepplr") # while developing the package, download from my fork
library(treepplr)
library(ape)
library(ggtree)
library(evolnets)
library(coda)
library(RevGadgets)
library(patchwork)
library(kdensity)
library(tidyverse)
```

## 1. Prepare data for analyses

Create symbiont tree, host tree, and interaction matrix.

```{r}
#| eval: false
#| label: simulate_data

tree <- rcoal(4, rooted = TRUE)
tree$tip.label <- paste0("S", 1:4)
height <- node.depth.edgelength(tree)[1]
tree$edge.length <- tree$edge.length * (2.0/height)

is.binary(tree)
is.ultrametric(tree)
is.rooted(tree)

host_tree <- rcoal(3, rooted = TRUE)
host_tree$tip.label <- paste0("H", 1:3)
plot(host_tree)
axisPhylo()

matrix <- matrix(data = c(2,0,2, 2,0,0, 2,2,0, 0,2,0), nrow = 4, ncol = 3, byrow = TRUE)
rownames(matrix) <- tree$tip.label
colnames(matrix) <- host_tree$tip.label

# write files
write.csv(matrix, "data/matrix.csv", row.names = TRUE)
write.nexus.data(matrix, "data/matrix.nex", format = "standard")
write.tree(host_tree, "data/host_tree.tre")
write.tree(tree, "data/tree.tre")

# add stem branch to symbiont tree
tree_string <- readLines("data/tree.tre")
tree_tiny_stem_string <- sub(");$", "):0.01;", tree_string)
tree_long_stem_string <- sub(");$", "):2.0;", tree_string)

writeLines(tree_tiny_stem_string, "data/tree_tiny_stem.tre")
writeLines(tree_long_stem_string, "data/tree_long_stem.tre")
```

```{r}
#| echo: false
#| label: read_data

tree <- read.tree("data/tree.tre")
host_tree <- read.tree("data/host_tree.tre")
matrix <- read.csv("data/matrix.csv", row.names = 1) %>% as.matrix()
```

Because RevBayes has it's own way to label tree nodes, we'll use the symbiont tree processed by RevBayes to keep track of node labels and be able to match them later. Then we can plot all data together.

```{r}
#| label: plot_data
#| fig-width: 7
#| fig-height: 6

symbiont_tree <- read_tree_from_revbayes("data/tree_tiny_stem_Rev.tre")

plot_data <- plot_matrix_phylo(matrix, at_nodes = NULL, symbiont_tree, host_tree, find_modules = FALSE)
plot_data[[1]] <- plot_data[[1]] + geom_nodelab(size = 3, hjust = -0.1)
plot_data

```


## 2. Run treeppl

```{r}
#| label: run_treeppl
#| eval: false

# symbiont_tree: TreeLabeled, ntips: Int, nhosts: Int,
# interactions: Int[], host_distances: Real[],
# dMean: Real, tune: Real

ntips <- Ntip(symbiont_tree)
nhosts <- Ntip(host_tree)
interactions <- matrix
host_distances <- cophenetic.phylo(host_tree)
dMean <- sum(host_distances) / factorial(nhosts)
tune <- 0.9

# I can't write a tree object to the json format that treeppl understands yet,
# so I'll run passing the json input file directly

out_treeppl <- run_treeppl(source = "/Users/mari/repos/treepplxrevbayes/code/treeppl/host_repertoire.tppl", 
                           data_path = "./data/test4x3_input.json",
                           method = "smc-apf",  
                           samples = 1000, chains = 100, sample_freq = 1000,
                           out_path = "./output/tiny_stem/treeppl/out_r_11.json")

```


## 3. Output analysis

### 3.1 RevBayes

```{r}
#| label: convergence_rb

long_rb_log1 <- read.table("output/long_stem/revbayes/long.1.test4x3.log", header = TRUE)
long_rb_log2 <- read.table("output/long_stem/revbayes/long.2.test4x3.log", header = TRUE)

# take only columns of interest
long_rb_log1 <- long_rb_log1[,c(1,5:6,8:11)] %>% filter(Iteration > 100000 & Iteration <= 700000)
long_rb_log2 <- long_rb_log2[,c(1,5:6,8:11)] %>% filter(Iteration > 100000 & Iteration <= 700000)

# give them better column names
colnames(long_rb_log1) <- colnames(long_rb_log2) <- c("iteration","clock","beta", "gain01", "loss10", "gain12", "loss21")

# convergence test
gelman.diag(mcmc.list(as.mcmc(long_rb_log1[,1:2]), as.mcmc(long_rb_log2[,1:2])))

effectiveSize(long_rb_log1)
effectiveSize(long_rb_log2)
```

Plot posterior distribution for each parameter

```{r}
#| label: param_dens_rb
#| fig-width: 5
#| fig-height: 7
#| eval: false
#| echo: false

trace_long_rb <- list(long_rb_log2)

params <- colnames(long_rb_log1)[-1]

dens_long_rb_list <- list()
for(i in seq_along(params)){
  dens_long_rb_list[i] <- plotTrace(trace_long_rb, vars = params[i])  
}
wrap_plots(dens_long_rb_list, ncol = 2)

```

Test if beta > 0

```{r}
#| label: bayes_factor

# Bayes factor

d_prior <- dexp(x=0, rate=1)

kd_beta <- kdensity(x = long_rb_log2$beta,
                    kernel='gamma',
                    support=c(0,Inf),
                    bw = 0.02)
max = kd_beta(0)

(BF <- d_prior/max)  # beta not different from zero
```


### 3.2 TreePPL

```{r}
files <- list.files(path = "./output/positive_debt", pattern = "*100k.json", full.names = TRUE)
files <- files[-3]

treeppl_logs <- data.frame()
for(i in seq_along(files)){
  log <- read_treeppl_output(files[i], sample_freq = 1000)
  treeppl_logs <- rbind(treeppl_logs, log)
}

treeppl_logs <- treeppl_logs %>% 
  mutate(run = 1:nrow(.)) %>% 
  mutate(posterior = exp(normConst - max(normConst)),
         sweep = dense_rank(normConst)) 

(zs <- sort(unique(treeppl_logs$normConst)))
var(zs)
```

Figure out how many subsamples to take

```{r}
treeppl_logs %>% 
  filter(norm_weight_within_sweep > 0.001) %>% 
  ggplot() +
  geom_density(aes(norm_weight_within_sweep, after_stat(scaled)), fill = "darkorange", col = NA, alpha = 0.7) +
  facet_wrap(vars(normConst)) +
  theme_bw()

treeppl_logs %>% 
  filter(norm_weight_within_sweep > 0.001) %>% 
  ggplot() +
  geom_col(aes(run, norm_weight_within_sweep), width = 0.1) +
  facet_wrap(vars(sweep), scales = "free")
  
treeppl_logs %>% 
  filter(sweep == 2) %>% 
  filter(norm_weight_within_sweep > 0.001) %>% 
  ggplot() +
  geom_col(aes(run, norm_weight_within_sweep), width = 0.1)

treeppl_logs %>% 
  filter(sweep == 1) %>% 
  pull(norm_weight_within_sweep) %>% 
  sort(decreasing = TRUE)

```


Plot parameter estimates from both RevBayes and TreePPL using the the normalized weights

```{r}
gg1 <- ggplot() +
  geom_density(aes(clock, after_stat(scaled)), fill = "steelblue", col = NA, alpha = 0.7, data = long_rb_log2) +
  geom_col(aes(mu, norm_weight_within_sweep), width = 0.01, alpha = 0.7, data = treeppl_logs) +
  #geom_density(aes(mu, weight = norm_weight_within_sweep, after_stat(scaled)), fill = "darkorange", col = NA, alpha = 0.7, data = treeppl_logs) +
  theme_bw()

gg2 <- ggplot() +
  geom_density(aes(beta, after_stat(scaled)), fill = "steelblue", col = NA, alpha = 0.7, data = long_rb_log2) +
  geom_col(aes(beta, norm_weight_within_sweep), width = 0.01, alpha = 0.7, data = filter(treeppl_logs, sweep == 1)) +
  #geom_density(aes(beta, weight = norm_weight_within_sweep, after_stat(scaled)), fill = "darkorange", col = NA, alpha = 0.7, data = treeppl_logs) +
  theme_bw()

gg3 <- ggplot() +
  geom_density(aes(gain01, after_stat(scaled)), fill = "steelblue", col = NA, alpha = 0.7, data = long_rb_log2) +
  geom_col(aes(lambda_1, norm_weight_within_sweep), width = 0.01, alpha = 0.7, data = treeppl_logs) +
  #geom_density(aes(lambda_1, weight = norm_weight_within_sweep, after_stat(scaled)), fill = "darkorange", col = NA, alpha = 0.7, data = treeppl_logs) +
  theme_bw()

gg4 <- ggplot() +
  geom_density(aes(loss10, after_stat(scaled)), fill = "steelblue", col = NA, alpha = 0.7, data = long_rb_log2) +
  geom_col(aes(lambda_2, norm_weight_within_sweep), width = 0.01, alpha = 0.7, data = treeppl_logs) +
  #geom_density(aes(lambda_2, weight = norm_weight_within_sweep, after_stat(scaled)), fill = "darkorange", col = NA, alpha = 0.7, data = treeppl_logs) +
  theme_bw()

gg5 <- ggplot() +
  geom_density(aes(gain12, after_stat(scaled)), fill = "steelblue", col = NA, alpha = 0.7, data = long_rb_log2) +
  geom_col(aes(lambda_3, norm_weight_within_sweep), width = 0.01, alpha = 0.7, data = treeppl_logs) +
  #geom_density(aes(lambda_3, weight = norm_weight_within_sweep, after_stat(scaled)), fill = "darkorange", col = NA, alpha = 0.7, data = treeppl_logs) +
  theme_bw()

gg6 <- ggplot() +
  geom_density(aes(loss21, after_stat(scaled)), fill = "steelblue", col = NA, alpha = 0.7, data = long_rb_log2) +
  geom_col(aes(lambda_4, norm_weight_within_sweep), width = 0.01, alpha = 0.7, data = treeppl_logs) +
  #geom_density(aes(lambda_4, weight = norm_weight_within_sweep, after_stat(scaled)), fill = "darkorange", col = NA, alpha = 0.7, data = treeppl_logs) +
  theme_bw()

gg1 + gg2 + gg3 + gg4 + gg5 + gg6 + plot_layout(ncol = 2, nrow = 3, byrow = TRUE)

```

Choose number of particles per sweep
```{r}
nsamples <- 3
subset <- treeppl_logs %>%
  group_by(sweep) %>% 
  slice_max(norm_weight_within_sweep, n = nsamples) %>% 
  mutate(wposterior = norm_weight_within_sweep * posterior)

```



Plot parameter estimates from both RevBayes and TreePPL using the normalizing constant

```{r}
g1 <- ggplot() +
  geom_density(aes(clock, after_stat(scaled)), fill = "steelblue", col = NA, alpha = 0.7, data = long_rb_log2) +
  geom_density(aes(mu, weight = posterior, after_stat(scaled)), fill = "darkorange", col = NA, alpha = 0.7, data = subset) +
  theme_bw()

g2 <- ggplot() +
  geom_density(aes(beta, after_stat(scaled)), fill = "steelblue", col = NA, alpha = 0.7, data = long_rb_log2) +
  geom_density(aes(beta, weight = posterior, after_stat(scaled)), fill = "darkorange", col = NA, alpha = 0.7, data = subset) +
  theme_bw()

g3 <- ggplot() +
  geom_density(aes(gain01, after_stat(scaled)), fill = "steelblue", col = NA, alpha = 0.7, data = long_rb_log2) +
  geom_density(aes(lambda_1, weight = posterior, after_stat(scaled)), fill = "darkorange", col = NA, alpha = 0.7, data = subset) +
  theme_bw()

g4 <- ggplot() +
  geom_density(aes(loss10, after_stat(scaled)), fill = "steelblue", col = NA, alpha = 0.7, data = long_rb_log2) +
  geom_density(aes(lambda_2, weight = posterior, after_stat(scaled)), fill = "darkorange", col = NA, alpha = 0.7, data = subset) +
  theme_bw()

g5 <- ggplot() +
  geom_density(aes(gain12, after_stat(scaled)), fill = "steelblue", col = NA, alpha = 0.7, data = long_rb_log2) +
  geom_density(aes(lambda_3, weight = posterior, after_stat(scaled)), fill = "darkorange", col = NA, alpha = 0.7, data = subset) +
  theme_bw()

g6 <- ggplot() +
  geom_density(aes(loss21, after_stat(scaled)), fill = "steelblue", col = NA, alpha = 0.7, data = long_rb_log2) +
  geom_density(aes(lambda_4, weight = posterior, after_stat(scaled)), fill = "darkorange", col = NA, alpha = 0.7, data = subset) +
  theme_bw()

g1 + g2 + g3 + g4 + g5 + g6 + plot_layout(ncol = 2, nrow = 3, byrow = TRUE)

```


Plot parameter estimates from both RevBayes and TreePPL using the normalizing constant * norm weight within sweep

```{r}
wg1 <- ggplot() +
  geom_density(aes(clock, after_stat(scaled)), fill = "steelblue", col = NA, alpha = 0.7, data = long_rb_log2) +
  geom_density(aes(mu, weight = wposterior, after_stat(scaled)), fill = "darkorange", col = NA, alpha = 0.7, data = subset) +
  theme_bw()

wg2 <- ggplot() +
  geom_density(aes(beta, after_stat(scaled)), fill = "steelblue", col = NA, alpha = 0.7, data = long_rb_log2) +
  geom_density(aes(beta, weight = wposterior, after_stat(scaled)), fill = "darkorange", col = NA, alpha = 0.7, data = subset) +
  theme_bw()

wg3 <- ggplot() +
  geom_density(aes(gain01, after_stat(scaled)), fill = "steelblue", col = NA, alpha = 0.7, data = long_rb_log2) +
  geom_density(aes(lambda_1, weight = wposterior, after_stat(scaled)), fill = "darkorange", col = NA, alpha = 0.7, data = subset) +
  theme_bw()

wg4 <- ggplot() +
  geom_density(aes(loss10, after_stat(scaled)), fill = "steelblue", col = NA, alpha = 0.7, data = long_rb_log2) +
  geom_density(aes(lambda_2, weight = wposterior, after_stat(scaled)), fill = "darkorange", col = NA, alpha = 0.7, data = subset) +
  theme_bw()

wg5 <- ggplot() +
  geom_density(aes(gain12, after_stat(scaled)), fill = "steelblue", col = NA, alpha = 0.7, data = long_rb_log2) +
  geom_density(aes(lambda_3, weight = wposterior, after_stat(scaled)), fill = "darkorange", col = NA, alpha = 0.7, data = subset) +
  theme_bw()

wg6 <- ggplot() +
  geom_density(aes(loss21, after_stat(scaled)), fill = "steelblue", col = NA, alpha = 0.7, data = long_rb_log2) +
  geom_density(aes(lambda_4, weight = wposterior, after_stat(scaled)), fill = "darkorange", col = NA, alpha = 0.7, data = subset) +
  theme_bw()

wg1 + wg2 + wg3 + wg4 + wg5 + wg6 + plot_layout(ncol = 2, nrow = 3, byrow = TRUE)

```


```{r}
ggplot() +
  geom_col(aes(lambda_1, posterior), width = 0.01, alpha = 0.7, data = subset, fill = "darkorange") +
  geom_col(aes(lambda_1, wposterior), width = 0.01, alpha = 0.7, data = subset, fill = "steelblue") +

  theme_bw()
```