Skip to content

Commit

Permalink
update model and prior root
Browse files Browse the repository at this point in the history
  • Loading branch information
sdgamboa committed Dec 5, 2023
1 parent b09c522 commit 105ab20
Show file tree
Hide file tree
Showing 2 changed files with 19 additions and 23 deletions.
20 changes: 11 additions & 9 deletions vignettes/articles/uncertainty_binary.Rmd
Original file line number Diff line number Diff line change
Expand Up @@ -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) |>
Expand All @@ -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)
```

Expand All @@ -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')
```
Expand All @@ -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')
```
Expand All @@ -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')
```
Expand Down
22 changes: 8 additions & 14 deletions vignettes/articles/uncertainty_multistate.Rmd
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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"))
Expand All @@ -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"))
Expand All @@ -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"))
Expand Down

0 comments on commit 105ab20

Please sign in to comment.