Uncertainty in tips with multistate(-intersection) annotations
Source:vignettes/articles/uncertainty_multistate.Rmd
uncertainty_multistate.Rmd
options(rstudio.viewer.autorefresh = FALSE)
library(ape)
library(phytools)
#> Loading required package: maps
library(dplyr)
#>
#> Attaching package: 'dplyr'
#> The following object is masked from 'package:ape':
#>
#> where
#> The following objects are masked from 'package:stats':
#>
#> filter, lag
#> The following objects are masked from 'package:base':
#>
#> intersect, setdiff, setequal, union
library(tidyr)Prepare data and tree
## 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
data <- primate.data
data <- data[tree$tip.label,]
original <- data |>
tibble::rownames_to_column(var = 'Taxa') |>
select(Taxa, Activity_pattern) |>
mutate(Presence = 1) |>
pivot_wider(
names_from = 'Activity_pattern', values_from = 'Presence',
values_fill = 0
) |>
arrange(Taxa) |>
tibble::column_to_rownames(var = 'Taxa') |>
as.matrix()
colnames(original) <- c(LETTERS[1:3])
myFun <- function(mat, per_uncertain = 0.1) {
n_row <- nrow(mat)
n <- round(n_row * per_uncertain)
rows <- sample(x = 1:nrow(mat), size = n, replace = FALSE)
mat[rows,] <- rep(1/ncol(mat), ncol(mat))
mat
}No uncertainty
fit <- fitMk(tree = tree, x = original,
model = myMod, pi = myPi,
lik.func = "pruning", logscale = TRUE)
ace <- ancr(fit, tips=TRUE)
plot(ace, args.plotTree = list(direction="upwards"))
title(main = '0% uncertain tips', line = -1)
Uncertainty about 12%
macaca <- grep('^Macaca_', rownames(original)) # A
galago <- grep('^Galago_', rownames(original)) # B
eulemur <- grep('^Eulemur', rownames(original)) # C
paste0(round((length(macaca) + length(galago) + length(eulemur)) / 90 * 100), '%')
#> [1] "12%"
m1 <- original
m1[] <- 1 / 3
m1[macaca, ] <- c(rep(1,3), rep(0, 6))
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 = myMod, pi = myPi,
lik.func = "pruning", logscale = TRUE)
ace1 <- ancr(fit1, tips=TRUE)
plot(ace1, args.plotTree = list(direction = "upwards"))
title(main = '12% uncertain tips', line = -1)
Let’s add some more taxa
cebus <- grep('^Cebus_', rownames(original)) # A
microcebus <- grep('^Microcebus', rownames(original)) # B
paste0(round((length(macaca) + length(galago) + length(eulemur) + length(cebus) + length(microcebus)) / 90 * 100), '%')
#> [1] "17%"
m2 <- original
m2[] <- 1 / 3
m2[macaca, ] <- c(rep(1,3), rep(0, 6))
m2[galago, ] <- c(rep(0,4), rep(1, 4), rep(0, 4))
m2[eulemur, ] <- c(rep(0,8), rep(1, 4))
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 = myMod, pi = myPi,
lik.func = "pruning", logscale = TRUE)
ace2 <- ancr(fit2, tips=TRUE)
plot(ace2, args.plotTree = list(direction = "upwards"))
title(main = '17% uncertain tips', line = -1)
Session information
sessioninfo::session_info()
#> ─ Session info ───────────────────────────────────────────────────────────────
#> setting value
#> version R version 4.4.1 (2024-06-14)
#> os Ubuntu 22.04.4 LTS
#> system x86_64, linux-gnu
#> ui X11
#> language en
#> collate en_US.UTF-8
#> ctype en_US.UTF-8
#> tz Etc/UTC
#> date 2024-11-21
#> pandoc 3.2 @ /usr/bin/ (via rmarkdown)
#>
#> ─ Packages ───────────────────────────────────────────────────────────────────
#> package * version date (UTC) lib source
#> ape * 5.8 2024-04-11 [1] RSPM (R 4.4.0)
#> bslib 0.8.0 2024-07-29 [1] RSPM (R 4.4.0)
#> cachem 1.1.0 2024-05-16 [1] RSPM (R 4.4.0)
#> cli 3.6.3 2024-06-21 [1] RSPM (R 4.4.0)
#> clusterGeneration 1.3.8 2023-08-16 [1] RSPM (R 4.4.0)
#> coda 0.19-4.1 2024-01-31 [1] RSPM (R 4.4.0)
#> codetools 0.2-20 2024-03-31 [2] CRAN (R 4.4.1)
#> combinat 0.0-8 2012-10-29 [1] RSPM (R 4.4.0)
#> DEoptim 2.2-8 2022-11-11 [1] RSPM (R 4.4.0)
#> desc 1.4.3 2023-12-10 [1] RSPM (R 4.4.0)
#> digest 0.6.37 2024-08-19 [1] RSPM (R 4.4.0)
#> doParallel 1.0.17 2022-02-07 [1] RSPM (R 4.4.0)
#> dplyr * 1.1.4 2023-11-17 [1] RSPM (R 4.4.0)
#> evaluate 1.0.1 2024-10-10 [1] RSPM (R 4.4.0)
#> expm 1.0-0 2024-08-19 [1] RSPM (R 4.4.0)
#> fansi 1.0.6 2023-12-08 [1] RSPM (R 4.4.0)
#> fastmap 1.2.0 2024-05-15 [1] RSPM (R 4.4.0)
#> fastmatch 1.1-4 2023-08-18 [1] RSPM (R 4.4.0)
#> foreach 1.5.2 2022-02-02 [1] RSPM (R 4.4.0)
#> fs 1.6.5 2024-10-30 [1] RSPM (R 4.4.0)
#> generics 0.1.3 2022-07-05 [1] RSPM (R 4.4.0)
#> glue 1.8.0 2024-09-30 [1] RSPM (R 4.4.0)
#> htmltools 0.5.8.1 2024-04-04 [1] RSPM (R 4.4.0)
#> htmlwidgets 1.6.4 2023-12-06 [1] RSPM (R 4.4.0)
#> igraph 2.1.1 2024-10-19 [1] RSPM (R 4.4.0)
#> iterators 1.0.14 2022-02-05 [1] RSPM (R 4.4.0)
#> jquerylib 0.1.4 2021-04-26 [1] RSPM (R 4.4.0)
#> jsonlite 1.8.9 2024-09-20 [1] RSPM (R 4.4.0)
#> knitr 1.49 2024-11-08 [1] RSPM (R 4.4.0)
#> lattice 0.22-6 2024-03-20 [2] CRAN (R 4.4.1)
#> lifecycle 1.0.4 2023-11-07 [1] RSPM (R 4.4.0)
#> magrittr 2.0.3 2022-03-30 [1] RSPM (R 4.4.0)
#> maps * 3.4.2.1 2024-11-10 [1] RSPM (R 4.4.0)
#> MASS 7.3-61 2024-06-13 [2] RSPM (R 4.4.0)
#> Matrix 1.7-0 2024-04-26 [2] CRAN (R 4.4.1)
#> mnormt 2.1.1 2022-09-26 [1] RSPM (R 4.4.0)
#> nlme 3.1-165 2024-06-06 [2] RSPM (R 4.4.0)
#> numDeriv 2016.8-1.1 2019-06-06 [1] RSPM (R 4.4.0)
#> optimParallel 1.0-2 2021-02-11 [1] RSPM (R 4.4.0)
#> phangorn 2.12.1 2024-09-17 [1] RSPM (R 4.4.0)
#> phytools * 2.3-0 2024-06-13 [1] RSPM (R 4.4.0)
#> pillar 1.9.0 2023-03-22 [1] RSPM (R 4.4.0)
#> pkgconfig 2.0.3 2019-09-22 [1] RSPM (R 4.4.0)
#> pkgdown 2.1.1 2024-09-17 [1] RSPM (R 4.4.0)
#> purrr 1.0.2 2023-08-10 [1] RSPM (R 4.4.0)
#> quadprog 1.5-8 2019-11-20 [1] RSPM (R 4.4.0)
#> R6 2.5.1 2021-08-19 [1] RSPM (R 4.4.0)
#> ragg 1.3.2 2024-05-15 [1] RSPM (R 4.4.0)
#> Rcpp 1.0.13-1 2024-11-02 [1] RSPM (R 4.4.0)
#> rlang 1.1.4 2024-06-04 [1] RSPM (R 4.4.0)
#> rmarkdown 2.29 2024-11-04 [1] RSPM (R 4.4.0)
#> sass 0.4.9 2024-03-15 [1] RSPM (R 4.4.0)
#> scatterplot3d 0.3-44 2023-05-05 [1] RSPM (R 4.4.0)
#> sessioninfo 1.2.2 2021-12-06 [1] RSPM (R 4.4.0)
#> systemfonts 1.1.0 2024-05-15 [1] RSPM (R 4.4.0)
#> textshaping 0.4.0 2024-05-24 [1] RSPM (R 4.4.0)
#> tibble 3.2.1 2023-03-20 [1] RSPM (R 4.4.0)
#> tidyr * 1.3.1 2024-01-24 [1] RSPM (R 4.4.0)
#> tidyselect 1.2.1 2024-03-11 [1] RSPM (R 4.4.0)
#> utf8 1.2.4 2023-10-22 [1] RSPM (R 4.4.0)
#> vctrs 0.6.5 2023-12-01 [1] RSPM (R 4.4.0)
#> withr 3.0.2 2024-10-28 [1] RSPM (R 4.4.0)
#> xfun 0.49 2024-10-31 [1] RSPM (R 4.4.0)
#> yaml 2.3.10 2024-07-26 [1] RSPM (R 4.4.0)
#>
#> [1] /usr/local/lib/R/site-library
#> [2] /usr/local/lib/R/library
#>
#> ──────────────────────────────────────────────────────────────────────────────