From 105ab209c9689b9b5e2a2f06e6d678d986ef6aff Mon Sep 17 00:00:00 2001 From: sdgamboa Date: Mon, 4 Dec 2023 21:52:55 -0500 Subject: [PATCH] update model and prior root --- vignettes/articles/uncertainty_binary.Rmd | 20 +++++++++-------- vignettes/articles/uncertainty_multistate.Rmd | 22 +++++++------------ 2 files changed, 19 insertions(+), 23 deletions(-) diff --git a/vignettes/articles/uncertainty_binary.Rmd b/vignettes/articles/uncertainty_binary.Rmd index b6968fb..b552bd5 100644 --- a/vignettes/articles/uncertainty_binary.Rmd +++ b/vignettes/articles/uncertainty_binary.Rmd @@ -35,6 +35,8 @@ data <- primate.data data <- data[tree$tip.label,] # rownames(data) <- paste0('taxon', 1:nrow(data)) # tree$tip.label <- paste0('taxon', 1:Ntip(tree)) +myMod <- 'ARD' # 'ARD' +myPi <- 'fitzjohn' original <- data |> tibble::rownames_to_column(var = 'Taxa') |> select(Taxa, Activity_pattern) |> @@ -60,10 +62,10 @@ head(sort(table(sub('^(\\w+)_.*$', '\\1', rownames(original))), decreasing = TRU ```{r} fit <- fitMk(tree = tree, x = original, - model = "ARD", pi = "fitzjohn", + model = myMod, pi = myPi, lik.func = "pruning", logscale = TRUE) -ace <- ancr(fit, tips=TRUE) -plot(ace, args.plotTree = list(direction="upwards")) +ace <- ancr(fit, tips = TRUE) +plot(ace, args.plotTree = list(direction = "upwards")) title(main = '0% uncertain tips', line = -1) ``` @@ -74,9 +76,9 @@ m1[which(grepl('^Macaca_', rownames(m1))), 1] <- 1 m1[which(!grepl('^Macaca_', rownames(m1))), 1] <- 0 m1[,2] <- ifelse(m1[,1] == 1, 0, 1) fit1 <- fitMk(tree = tree, x = m1, - model = "ARD", pi = "fitzjohn", + model = myMod, pi = myPi, lik.func = "pruning", logscale = TRUE) -ace1 <- ancr(fit1, tips=TRUE) +ace1 <- ancr(fit1, tips = TRUE) plot(ace1, args.plotTree = list(direction = "upwards")) title(main = 'Macaca TRUE', line = -1, sub = '0/1 T/F for unknowns') ``` @@ -87,9 +89,9 @@ m2[] <- 0.5 m2[which(grepl('^Macaca_', rownames(m2))), 1] <- 1 m2[,2] <- ifelse(m2[,1] == 1, 0, 0.5) fit2 <- fitMk(tree = tree, x = m2, - model = "ARD", pi = "fitzjohn", + model = myMod, pi = myPi, lik.func = "pruning", logscale = TRUE) -ace2 <- ancr(fit2, tips=TRUE) +ace2 <- ancr(fit2, tips = TRUE) plot(ace2, args.plotTree = list(direction = "upwards")) title(main = 'Macaca TRUE', line = -1, sub = '0.5 T/F for unknowns') ``` @@ -103,9 +105,9 @@ m3[which(grepl('^Macaca_', rownames(m3))), 2] <- 0 m3[which(grepl('^Galago_', rownames(m3))), 2] <- 1 m3[which(grepl('^Galago', rownames(m3))), 1] <- 0 fit3 <- fitMk(tree = tree, x = m3, - model = "ARD", pi = "fitzjohn", + model = myMod, pi = myPi, lik.func = "pruning", logscale = TRUE) -ace3 <- ancr(fit3, tips=TRUE) +ace3 <- ancr(fit3, tips = TRUE) plot(ace3, args.plotTree = list(direction = "upwards")) title(main = 'Macaca TRUE - Galago FALSE', line = -1, sub = '0.5 T/F for unknowns') ``` diff --git a/vignettes/articles/uncertainty_multistate.Rmd b/vignettes/articles/uncertainty_multistate.Rmd index 892774b..8ab20b6 100644 --- a/vignettes/articles/uncertainty_multistate.Rmd +++ b/vignettes/articles/uncertainty_multistate.Rmd @@ -20,6 +20,11 @@ library(tidyr) ## Prepare data and tree ```{r} +## ARD with fitzjohn +## ER with fitzjohn +## ER with estimated +myMod <- 'ARD' # ARD*, ER +myPi <- 'fitzjohn' # fitzjohn*, equal, estimated data('primate.tree') data('primate.data') tree <- primate.tree @@ -49,22 +54,11 @@ myFun <- function(mat, per_uncertain = 0.1) { ``` -## Some practical examples - - -```{r} - -``` - - - - - # No uncertainty ```{r} fit <- fitMk(tree = tree, x = original, - model = "ARD", pi = "fitzjohn", + model = myMod, pi = myPi, lik.func = "pruning", logscale = TRUE) ace <- ancr(fit, tips=TRUE) plot(ace, args.plotTree = list(direction="upwards")) @@ -90,7 +84,7 @@ m1[galago, ] <- c(rep(0,4), rep(1, 4), rep(0, 4)) m1[eulemur, ] <- c(rep(0,8), rep(1, 4)) fit1 <- fitMk(tree = tree, x = m1, - model = "ARD", pi = "fitzjohn", + model = myMod, pi = myPi, lik.func = "pruning", logscale = TRUE) ace1 <- ancr(fit1, tips=TRUE) plot(ace1, args.plotTree = list(direction = "upwards")) @@ -117,7 +111,7 @@ m2[cebus, ] <- c(rep(1,2), rep(0, 4)) m2[microcebus, ] <- c(rep(0,2), rep(1, 2), rep(0, 2)) fit2 <- fitMk(tree = tree, x = m2, - model = "ARD", pi = "fitzjohn", + model = myMod, pi = myPi, lik.func = "pruning", logscale = TRUE) ace2 <- ancr(fit2, tips=TRUE) plot(ace2, args.plotTree = list(direction = "upwards"))