add interna_nodes.Rmd vignette/artticle

vignettes/articles/internal_nodes.Rmd

@@ -0,0 +1,154 @@
+---
+title: "internal_nodes"
+---
+
+```{r, include = FALSE}
+knitr::opts_chunk$set(
+  collapse = TRUE,
+  comment = "#>"
+)
+```
+
+```{r setup, message=FALSE}
+library(phytools)
+library(dplyr)
+library(tidyr)
+library(tibble)
+```
+
+## Data
+
+```{r}
+data('primate.tree')
+data('primate.data')
+tre <- primate.tree
+dat <- primate.data
+states_original <- dat |> 
+    rownames_to_column(var = 'tip_label') |> 
+    select(tip_label, Activity_pattern) |> 
+    mutate(state = 1) |> 
+    pivot_wider(
+        names_from = 'Activity_pattern', values_from = 'state', values_fill = 0
+    ) |> 
+    column_to_rownames(var = 'tip_label') |> 
+    as.matrix() |> 
+    magrittr::set_colnames(LETTERS[1:3])
+states_original <- states_original[tre$tip.label,]
+
+
+macaca <- grep('^Macaca_', rownames(states_original)) # ?
+galago <- grep('^Galago_', rownames(states_original)) # B
+eulemur <- grep('^Eulemur', rownames(states_original)) # C
+
+states <- states_original
+states[] <- 1 / ncol(states)
+# m1[macaca, ] <- c(rep(1,3), rep(0, 6))
+states[galago, ] <- c(rep(0,4), rep(1, 4), rep(0, 4))
+states[eulemur, ] <- c(rep(0,8), rep(1, 4))
+myMod <- 'ARD'
+myPi <- 'fitzjohn'
+```
+
+
+## Annotate one genus with an additional tip of zero length
+
+This is the recommended way of adding an annotations to an internal node.
+
+```{r}
+macaca <- grep('Macaca', rownames(states))
+macaca_mrc <- findMRCA(tree = tre, tips = rownames(states)[macaca])
+t1 <- bind.tip(tre, tip.label = 'Macaca_gn1', edge.length = 0, where = macaca_mrc)
+m1 <- states_original
+# m1[] <- 1 / 3
+m1 <- rbind(m1, matrix(c(A = 0, B = 0, C = 1), nrow = 1, byrow = TRUE))
+rownames(m1)[nrow(m1)] <- 'Macaca_gn1'
+m1[macaca,] <- rep(1 / ncol(m1), length(macaca) * 3)
+
+fit1 <- fitMk(
+    tree = t1, x = m1, model = myMod, pi = myPi, lik.func = "pruning",
+    logscale = TRUE
+)
+ace1 <- ancr(fit1, tips=TRUE)
+plot(ace1, args.plotTree = list(direction = "upwards"))
+```
+
+However, in some cases, one might need to add two or more tips carrying
+conclicting annotations. If both are of length zero, this throws an error:
+
+```{r}
+macaca_mrc2 <- findMRCA(t1, tips = rownames(states)[macaca])
+t2 <- bind.tip(t1, tip.label = 'Macaca_gn2', edge.length = 0, where = macaca_mrc2)
+m2 <- states
+m2 <- rbind(m2, matrix(c(A = 0, B = 0, C = 1), nrow = 1, byrow = TRUE))
+m2 <- rbind(m2, matrix(c(A = 0, B = 1, C = 0), nrow = 1, byrow = TRUE))
+rownames(m2)[nrow(m2) - 1] <- 'Macaca_gn1'
+rownames(m2)[nrow(m2)] <- 'Macaca_gn2'
+m2[grep('^Macaca_[^g][^n][^\\d].*', rownames(m2)),] <- rep(1 / ncol(m2), length(macaca) * 3)
+fit2 <- tryCatch(
+    error = function(e) e, 
+    {
+        fitMk(
+            tree = t2, x = m2, model = myMod, pi = myPi,
+            lik.func = "pruning", logscale = TRUE
+        )
+    }
+)
+print(fit2$message)
+```
+
+## Annotate one genus with an additional tip of length 1e-06
+
+A solultion to the issue described above could be adding small quantity to the
+branch length of the additional tips instead of a hard zero, e.g., 1e-06. At
+least with the example below, the output seems to be the same (compare with
+the first very firt case in this document).
+
+```{r}
+macaca <- grep('Macaca', rownames(states))
+macaca_mrc <- findMRCA(tree = tre, tips = rownames(states)[macaca])
+t1 <- bind.tip(tre, tip.label = 'Macaca_gn1', edge.length = 1e-06, where = macaca_mrc)
+m1 <- states_original
+# m1[] <- 1 / 3
+m1 <- rbind(m1, matrix(c(A = 0, B = 0, C = 1), nrow = 1, byrow = TRUE))
+rownames(m1)[nrow(m1)] <- 'Macaca_gn1'
+m1[macaca,] <- rep(1 / ncol(m1), length(macaca) * 3)
+
+fit1 <- fitMk(
+    tree = t1, x = m1, model = myMod, pi = myPi, lik.func = "pruning",
+    logscale = TRUE
+)
+ace1 <- ancr(fit1, tips=TRUE)
+plot(ace1, args.plotTree = list(direction = "upwards"))
+```
+
+## Addint two tips with branch length of 1e-06
+
+Now let's see if this has the expected result if I have two annotations.
+For this, I'll make all annotations uncertain, except for the tips of the genus 
+Eulermur (C/red) and Galago (B/white). I'll add two conflicting annotations
+to the genus Macaca using two tips with branch length of 1e-06.
+This seems to have given the expected result: half and half.
+
+```{r}
+t3 <- bind.tip(tre, tip.label = 'Macaca_gn1', edge.length = 1e-06, where = macaca_mrc)
+macaca_mrc3 <- findMRCA(t3, tips = rownames(states)[macaca])
+t4 <- bind.tip(t3, tip.label = 'Macaca_gn2', edge.length = 1e-06, where = macaca_mrc3)
+m3 <- states
+m3 <- rbind(m3, matrix(c(A = 1, B = 0, C = 0), nrow = 1, byrow = TRUE))
+m3 <- rbind(m3, matrix(c(A = 0, B = 1, C = 0), nrow = 1, byrow = TRUE))
+rownames(m3)[nrow(m3) - 1] <- 'Macaca_gn1'
+rownames(m3)[nrow(m3)] <- 'Macaca_gn2'
+m3[grep('^Macaca_[^g][^n][^\\d].*', rownames(m3)),] <- rep(1 / ncol(m3), length(macaca) * 3)
+fit3 <- fitMk(
+    tree = t4, x = m3, model = myMod, pi = myPi, lik.func = "pruning",
+    logscale = TRUE
+)
+ace3 <- ancr(fit3, tips=TRUE)
+plot(ace3, args.plotTree = list(direction = "upwards"))
+```
+
+# Session information
+
+```{r}
+sessioninfo::session_info()
+```