joshwlambert / daisieprep Goto Github PK

View Code? Open in Web Editor NEW

6.0 6.0 3.0 14.09 MB

Extracts phylogenetic island community data from phylogenetic trees

Home Page: https://joshwlambert.github.io/DAISIEprep

License: GNU General Public License v3.0

R 97.03% TeX 2.97%

data-science island-biogeography phylogenetics r

daisieprep's Introduction

Hi 👋

I'm Josh, an evolutionary biologist 🧬 🐸, turned research software engineer 💻.

I'm currently working at the The London School of Hygiene & Tropical Medicine (LSHTM) on the Epiverse initiative, and finishing my PhD at the University of Groningen. My PhD is on the topic of developing novel phylogenetic models to understand macroevolutionary patterns on islands ⛰️.

🔭 I’m currently working on: Epiverse by day, macroevolution and island biogeography by night.
🌱 I’m currently learning: most things, but mainly epidemiology of infectious diseases.
👯 I’m looking to collaborate on: interesting research questions at the intersection of evolutionary biology and epidemiology.
🤔 I’m looking for help with: not sure ... but fun collaborators are always welcome 😃
💬 Ask me about: anything!
📫 How to reach me: [email protected]
😄 Pronouns: he/him
⚡ Fun fact: I'm a fan of all things Monocle

daisieprep's People

Stargazers

Watchers

Forkers

neves-p theopannetier bincente

daisieprep's Issues

Improve documentation for `add_island_colonist()`

Comments helpfully provided by @rsetienne:

For status, please indicate what the options are.
For missing species, please indicate that for a clade with only missing species this should be one less than the number of missing species.
For clade_type, please indicate what this is. Is this for a type_2 analysis, or something else?

Add error message when add species to non-existing clade

When using the function add_missing_species() one can specify which clade to assign it to using: e.g. species_to_add_to = "Plant_e". This works fine. However, I noticed that if I add to as species that is not included in the island_tbl species_to_add_to = "Plant_that_is_not_in_the_islandtbl" the code runs and nothing happens. I think in this case it should give the following error message:

"You are adding species to a clade that is already sampled in the data. However, in species_to_add_to you are giving the name of a species that does not occur in the data, so the clade cannot be located."

Consider switching endemicity to state encoding

As of now, following the Tutorial and the MuSSE example, we have the following encoding:

not_present = 1
nonendemic = 2
endemic = 3

This is arbitrary, works well. I have hard-encoded the convention into all_endemicity_status().

However, GeoSSE requires a different encoding:

Area A only (mainland, not_present) = 1
Area B only (island, endemic) = 2
Area AB (widespread state,nonendemic) = 3

The output of BioGeoBEARS also appears to follow this convention. I'll switch the current convention to this one to avoid the coexistence of two encodings and the risk of confusion.

`extract_island_species` also extracts a mainland species?

Running code from the Tutorial vignette:

library(DAISIEprep)
library(ape)
library(phylobase)
#> 
#> Attaching package: 'phylobase'
#> The following object is masked from 'package:ape':
#> 
#>     edges
library(ggtree)
#> ggtree v3.5.1.900 For help: https://yulab-smu.top/treedata-book/
#> 
#> If you use the ggtree package suite in published research, please cite
#> the appropriate paper(s):
#> 
#> Guangchuang Yu, David Smith, Huachen Zhu, Yi Guan, Tommy Tsan-Yuk Lam.
#> ggtree: an R package for visualization and annotation of phylogenetic
#> trees with their covariates and other associated data. Methods in
#> Ecology and Evolution. 2017, 8(1):28-36. doi:10.1111/2041-210X.12628
#> 
#> G Yu. Data Integration, Manipulation and Visualization of Phylogenetic
#> Trees (1st ed.). Chapman and Hall/CRC. 2022. ISBN: 9781032233574
#> 
#> S Xu, Z Dai, P Guo, X Fu, S Liu, L Zhou, W Tang, T Feng, M Chen, L
#> Zhan, T Wu, E Hu, Y Jiang, X Bo, G Yu. ggtreeExtra: Compact
#> visualization of richly annotated phylogenetic data. Molecular Biology
#> and Evolution. 2021, 38(9):4039-4042. doi: 10.1093/molbev/msab166
#> 
#> Attaching package: 'ggtree'
#> The following object is masked from 'package:phylobase':
#> 
#>     MRCA
#> The following object is masked from 'package:ape':
#> 
#>     rotate
library(ggimage)
#> Loading required package: ggplot2
library(castor)
#> Loading required package: Rcpp
set.seed(
  1,
  kind = "Mersenne-Twister",
  normal.kind = "Inversion",
  sample.kind = "Rejection"
)
phylo <- ape::rcoal(10)
phylo$tip.label <- c("bird_a", "bird_b", "bird_c", "bird_d", "bird_e", "bird_f",
                     "bird_g", "bird_h", "bird_i", "bird_j")
phylo <- phylobase::phylo4(phylo)
endemicity_status <- sample(
  x = c("not_present", "endemic", "nonendemic"),
  size = length(phylobase::tipLabels(phylo)),
  replace = TRUE,
  prob = c(0.6, 0.2, 0.2)
)
phylod <- phylobase::phylo4d(phylo, as.data.frame(endemicity_status))
island_tbl <- island_tbl()
tip_states <- c()
endemicity_status <- phylobase::tipData(phylod)$endemicity_status
for (i in seq_along(endemicity_status)) {
  if (grepl(pattern = "^not_present$", x = endemicity_status[i])) {
    tip_states[i] <- 1
  } else if (grepl(pattern = "^nonendemic$", x = endemicity_status[i])) {
    tip_states[i] <- 2
  } else if (grepl(pattern = "^endemic$", x = endemicity_status[i])) {
    tip_states[i] <- 3
  }
}
phylo <- as(phylo, "phylo")
asr <- castor::asr_max_parsimony(
  tree = phylo,
  tip_states = tip_states,
  transition_costs = "sequential"
)
if (ncol(asr$ancestral_likelihoods) == 2) {
  colnames(asr$ancestral_likelihoods) <- c("not_present", "nonendemic")
} else if (ncol(asr$ancestral_likelihoods) == 3) {
  colnames(asr$ancestral_likelihoods) <-
    c("not_present", "nonendemic", "endemic")
}
node_states <- max.col(asr$ancestral_likelihoods, ties.method = "last")
node_states <- gsub(
  pattern = "1", replacement = "not_present", x = node_states
)
node_states <- gsub(
  pattern = "2", replacement = "nonendemic", x = node_states
)
node_states <- gsub(
  pattern = "3", replacement = "endemic", x = node_states
)
node_data <- data.frame(
  island_status = node_states,
  endemic_prob = asr$ancestral_likelihoods[, "endemic"],
  nonendemic_prob = asr$ancestral_likelihoods[, "nonendemic"],
  not_present_prob = asr$ancestral_likelihoods[, "not_present"],
  row.names = phylobase::nodeId(phylod, "internal")
)
phylod <- phylobase::phylo4d(
  phylo,
  tip.data = as.data.frame(endemicity_status),
  node.data = node_data
)
phylod <- add_asr_node_states(phylod = phylod, asr_method = "parsimony")
island_tbl <- extract_island_species(
  phylod = phylod,
  extraction_method = "asr",
  include_not_present = FALSE
)
island_tbl
#> Class:  Island_tbl 
#>   clade_name  status missing_species branching_times min_age      species
#> 1     bird_g endemic               0    0.764855....      NA bird_g, ....
island_tbl@island_tbl$species
#> [[1]]
#> [1] "bird_g" "bird_h" "bird_i"
island_tbl@island_tbl$branching_times
#> [[1]]
#> [1] 0.7648553 0.3800300
plot_phylod(phylod = phylod)

Shouldn't bird_h be left out by extract_island_species() when include_not_present = FALSE ?

^{Created on 2022-06-29 by the reprex package (v2.0.1)}

Duplicated extracted species from `extract_island_species()`

When a species, for this example let's call this species "bird_a", is extracted using the asr algorithm and if first stored as a colonist, it can then be duplicated in the island_tbl if another clade has an ancestral node on the island that also includes "bird_a" as one of it's descendants. It is not excluded in the same way a descendant that is not present on the island would be.

Here is a reprex to show the bug:

library(DAISIEprep)
set.seed(
  1,
  kind = "Mersenne-Twister",
  normal.kind = "Inversion",
  sample.kind = "Rejection"
)
phylo <- ape::rcoal(6)
phylo$tip.label <- c("bird_a", "bird_b", "bird_c", "bird_d", "bird_e",
                     "bird_f")
phylo <- add_outgroup(phylo = phylo)
phylo <- phylobase::phylo4(phylo)

endemicity_status <- c(
  "not_present",
  "endemic",
  "not_present",
  "nonendemic",
  rep("endemic", 2),
  "not_present"
)
phylod <- phylobase::phylo4d(phylo, as.data.frame(endemicity_status))

plot_phylod(phylod)

phylod <- add_asr_node_states(
  phylod = phylod, 
  asr_method = "mk",
  tie_preference = "mainland",
  rate_model = matrix(
    data = c(
      0, 1, 0,
      2, 0, 3,
      0, 0, 0
    ),
    nrow = 3,
    byrow = TRUE
  )
)

plot_phylod(phylod)

# augment node state of bird a ancestor to not present to reproduce bug
phylod@data$island_status[11] <- "not_present"

plot_phylod(phylod)

it <- extract_island_species(phylod, extraction_method = "asr")
it
#> Class:  Island_tbl 
#>   clade_name  status missing_species   col_time col_max_age branching_times
#> 1     bird_a endemic               0 0.05034546       FALSE              NA
#> 2     bird_c endemic               0 0.67546123       FALSE    0.239392....
#>   min_age      species clade_type
#> 1      NA       bird_a          1
#> 2      NA bird_a, ....          1
it@island_tbl$species
#> [[1]]
#> [1] "bird_a"
#> 
#> [[2]]
#> [1] "bird_a" "bird_c" "bird_d" "bird_e"

^{Created on 2024-05-06 with reprex v2.1.0}

Add Theo as a collaborator on the repo

Thank you :)

(I'm sure you know this but just in case -> Inviting collaborators to a personal repo)

Incorrect install instructions

Thanks for spotting @luislvalente!

Enable more ASR options from {castor}

Currently add_asr_node_states() only allows "parsimony" or "mk" as options to call either castor::asr_max_parsimony() or castor::asr_mk_model(), respectively. Adding a rate_model argument to add_asr_node_states() will allow more MK models to be selecting rather than fixing to the default equal-rates models.

Exported functions need better examples

Some exported functions do not have examples, others are outdated.

Read ASR doc

Scale bar in phylogenies generated with plot_phylod()

The scale bar in time calibrated phylogenies using plot_phylod() is oddly formatted. It would be nice if the present was time 0 and that the bar would end aligned with the tips.

Draft ASR vignette diversitree example

MuSSE example
MuSSE with constrains example
GeoSSE example

Draft ASR vignette corHMM example

`plot_phylod` is no longer working

I re-installed the package (from devtools::install(".")), and since reinstalling plot_phylod() keeps outputting this uninformative error:

Error in if (attr(treeview$data, "revts.done")) {: argument is of length zero

I haven't changed anything to the function, so I strongly suspect something broke with the dependencies. I updated a number of them during reinstall, but didn't note which ones.

@joshwlambert do you spot anything from the updated packages in the session info below / have you encountered the same issue?

library(DAISIEprep)
library(ape)
library(phylobase)
#> 
#> Attaching package: 'phylobase'
#> The following object is masked from 'package:ape':
#> 
#>     edges
library(ggtree)
#> ggtree v3.5.1.902 For help: https://yulab-smu.top/treedata-book/
#> 
#> If you use the ggtree package suite in published research, please cite
#> the appropriate paper(s):
#> 
#> Guangchuang Yu, David Smith, Huachen Zhu, Yi Guan, Tommy Tsan-Yuk Lam.
#> ggtree: an R package for visualization and annotation of phylogenetic
#> trees with their covariates and other associated data. Methods in
#> Ecology and Evolution. 2017, 8(1):28-36. doi:10.1111/2041-210X.12628
#> 
#> S Xu, Z Dai, P Guo, X Fu, S Liu, L Zhou, W Tang, T Feng, M Chen, L
#> Zhan, T Wu, E Hu, Y Jiang, X Bo, G Yu. ggtreeExtra: Compact
#> visualization of richly annotated phylogenetic data. Molecular Biology
#> and Evolution. 2021, 38(9):4039-4042. doi: 10.1093/molbev/msab166
#> 
#> Guangchuang Yu. Using ggtree to visualize data on tree-like structures.
#> Current Protocols in Bioinformatics. 2020, 69:e96. doi:10.1002/cpbi.96
#> 
#> Attaching package: 'ggtree'
#> The following object is masked from 'package:phylobase':
#> 
#>     MRCA
#> The following object is masked from 'package:ape':
#> 
#>     rotate
library(ggimage)
#> Loading required package: ggplot2
library(castor)
#> Loading required package: Rcpp
set.seed(
  1,
  kind = "Mersenne-Twister",
  normal.kind = "Inversion",
  sample.kind = "Rejection"
)
phylo <- ape::rcoal(10)
phylo$tip.label <- c("bird_a", "bird_b", "bird_c", "bird_d", "bird_e", "bird_f",
                     "bird_g", "bird_h", "bird_i", "bird_j")
phylo <- phylobase::phylo4(phylo)
endemicity_status <- sample(
  x = c("not_present", "endemic", "nonendemic"),
  size = length(phylobase::tipLabels(phylo)),
  replace = TRUE,
  prob = c(0.6, 0.2, 0.2)
)
phylod <- phylobase::phylo4d(phylo, as.data.frame(endemicity_status))
island_tbl <- island_tbl()
tip_states <- c()
endemicity_status <- phylobase::tipData(phylod)$endemicity_status
for (i in seq_along(endemicity_status)) {
  if (grepl(pattern = "^not_present$", x = endemicity_status[i])) {
    tip_states[i] <- 1
  } else if (grepl(pattern = "^nonendemic$", x = endemicity_status[i])) {
    tip_states[i] <- 2
  } else if (grepl(pattern = "^endemic$", x = endemicity_status[i])) {
    tip_states[i] <- 3
  }
}
phylo <- as(phylo, "phylo")
asr <- castor::asr_max_parsimony(
  tree = phylo,
  tip_states = tip_states,
  transition_costs = "sequential"
)
if (ncol(asr$ancestral_likelihoods) == 2) {
  colnames(asr$ancestral_likelihoods) <- c("not_present", "nonendemic")
} else if (ncol(asr$ancestral_likelihoods) == 3) {
  colnames(asr$ancestral_likelihoods) <-
    c("not_present", "nonendemic", "endemic")
}
node_states <- max.col(asr$ancestral_likelihoods, ties.method = "last")
node_states <- gsub(
  pattern = "1", replacement = "not_present", x = node_states
)
node_states <- gsub(
  pattern = "2", replacement = "nonendemic", x = node_states
)
node_states <- gsub(
  pattern = "3", replacement = "endemic", x = node_states
)
node_data <- data.frame(
  island_status = node_states,
  endemic_prob = asr$ancestral_likelihoods[, "endemic"],
  nonendemic_prob = asr$ancestral_likelihoods[, "nonendemic"],
  not_present_prob = asr$ancestral_likelihoods[, "not_present"],
  row.names = phylobase::nodeId(phylod, "internal")
)
phylod <- phylobase::phylo4d(
  phylo,
  tip.data = as.data.frame(endemicity_status),
  node.data = node_data
)
phylod <- add_asr_node_states(phylod = phylod, asr_method = "parsimony")
island_tbl <- extract_island_species(
  phylod = phylod,
  extraction_method = "asr",
  include_not_present = FALSE
)
island_tbl
#> Class:  Island_tbl 
#>   clade_name  status missing_species  col_time col_max_age branching_times
#> 1     bird_g endemic               0 0.7648553       FALSE    0.380034....
#>   min_age      species clade_type
#> 1      NA bird_g, ....          1
island_tbl@island_tbl$species
#> [[1]]
#> [1] "bird_g" "bird_i"
island_tbl@island_tbl$branching_times
#> [[1]]
#> [1] 0.3800341
plot_phylod(phylod = phylod)
#> Error in if (attr(treeview$data, "revts.done")) {: argument is of length zero

^{Created on 2022-07-25 by the reprex package (v2.0.1)}

Session info

sessioninfo::session_info()
#> ─ Session info ───────────────────────────────────────────────────────────────
#>  setting  value
#>  version  R version 4.1.2 (2021-11-01)
#>  os       macOS Catalina 10.15.6
#>  system   x86_64, darwin17.0
#>  ui       X11
#>  language (EN)
#>  collate  en_US.UTF-8
#>  ctype    en_US.UTF-8
#>  tz       Europe/Paris
#>  date     2022-07-25
#>  pandoc   2.17.1.1 @ /Applications/RStudio.app/Contents/MacOS/quarto/bin/ (via rmarkdown)
#> 
#> ─ Packages ───────────────────────────────────────────────────────────────────
#>  package      * version   date (UTC) lib source
#>  ade4           1.7-19    2022-04-19 [1] CRAN (R 4.1.2)
#>  ape          * 5.6-2     2022-03-02 [1] CRAN (R 4.1.2)
#>  aplot          0.1.6     2022-06-03 [1] CRAN (R 4.1.2)
#>  assertthat     0.2.1     2019-03-21 [1] CRAN (R 4.1.0)
#>  castor       * 1.7.2     2022-01-08 [1] CRAN (R 4.1.2)
#>  cli            3.3.0     2022-04-25 [1] CRAN (R 4.1.2)
#>  colorspace     2.0-3     2022-02-21 [1] CRAN (R 4.1.2)
#>  crayon         1.5.1     2022-03-26 [1] CRAN (R 4.1.2)
#>  DAISIEprep   * 0.1.0     2022-07-25 [1] local
#>  DBI            1.1.2     2021-12-20 [1] CRAN (R 4.1.0)
#>  digest         0.6.29    2021-12-01 [1] CRAN (R 4.1.0)
#>  dplyr          1.0.9     2022-04-28 [1] CRAN (R 4.1.2)
#>  ellipsis       0.3.2     2021-04-29 [1] CRAN (R 4.1.0)
#>  evaluate       0.15      2022-02-18 [1] CRAN (R 4.1.2)
#>  fansi          1.0.3     2022-03-24 [1] CRAN (R 4.1.2)
#>  fastmap        1.1.0     2021-01-25 [1] CRAN (R 4.1.0)
#>  fs             1.5.2     2021-12-08 [1] CRAN (R 4.1.0)
#>  generics       0.1.3     2022-07-05 [1] CRAN (R 4.1.2)
#>  ggfun          0.0.6     2022-04-01 [1] CRAN (R 4.1.2)
#>  ggimage      * 0.3.1     2022-04-25 [1] CRAN (R 4.1.2)
#>  ggplot2      * 3.3.6     2022-05-03 [1] CRAN (R 4.1.2)
#>  ggplotify      0.1.0     2021-09-02 [1] CRAN (R 4.1.0)
#>  ggtree       * 3.5.1.902 2022-07-25 [1] Github (YuLab-SMU/ggtree@19d342b)
#>  glue           1.6.2     2022-02-24 [1] CRAN (R 4.1.2)
#>  gridGraphics   0.5-1     2020-12-13 [1] CRAN (R 4.1.0)
#>  gtable         0.3.0     2019-03-25 [1] CRAN (R 4.1.0)
#>  highr          0.9       2021-04-16 [1] CRAN (R 4.1.0)
#>  hms            1.1.1     2021-09-26 [1] CRAN (R 4.1.0)
#>  htmltools      0.5.2     2021-08-25 [1] CRAN (R 4.1.0)
#>  httr           1.4.3     2022-05-04 [1] CRAN (R 4.1.2)
#>  jsonlite       1.8.0     2022-02-22 [1] CRAN (R 4.1.2)
#>  knitr          1.39      2022-04-26 [1] CRAN (R 4.1.2)
#>  lattice        0.20-45   2021-09-22 [1] CRAN (R 4.1.2)
#>  lazyeval       0.2.2     2019-03-15 [1] CRAN (R 4.1.0)
#>  lifecycle      1.0.1     2021-09-24 [1] CRAN (R 4.1.0)
#>  magick         2.7.3     2021-08-18 [1] CRAN (R 4.1.0)
#>  magrittr       2.0.3     2022-03-30 [1] CRAN (R 4.1.2)
#>  MASS           7.3-54    2021-05-03 [1] CRAN (R 4.1.2)
#>  Matrix         1.4-0     2021-12-08 [1] CRAN (R 4.1.0)
#>  munsell        0.5.0     2018-06-12 [1] CRAN (R 4.1.0)
#>  naturalsort    0.1.3     2016-08-30 [1] CRAN (R 4.1.0)
#>  nlme           3.1-153   2021-09-07 [1] CRAN (R 4.1.2)
#>  patchwork      1.1.1     2020-12-17 [1] CRAN (R 4.1.0)
#>  phylobase    * 0.8.10    2020-03-01 [1] CRAN (R 4.1.0)
#>  pillar         1.8.0     2022-07-18 [1] CRAN (R 4.1.2)
#>  pkgconfig      2.0.3     2019-09-22 [1] CRAN (R 4.1.0)
#>  plyr           1.8.7     2022-03-24 [1] CRAN (R 4.1.2)
#>  prettyunits    1.1.1     2020-01-24 [1] CRAN (R 4.1.0)
#>  progress       1.2.2     2019-05-16 [1] CRAN (R 4.1.0)
#>  purrr          0.3.4     2020-04-17 [1] CRAN (R 4.1.0)
#>  R.cache        0.15.0    2021-04-30 [1] CRAN (R 4.1.0)
#>  R.methodsS3    1.8.2     2022-06-13 [1] CRAN (R 4.1.2)
#>  R.oo           1.25.0    2022-06-12 [1] CRAN (R 4.1.2)
#>  R.utils        2.11.0    2021-09-26 [1] CRAN (R 4.1.0)
#>  R6             2.5.1     2021-08-19 [1] CRAN (R 4.1.0)
#>  Rcpp         * 1.0.9     2022-07-08 [1] CRAN (R 4.1.2)
#>  reprex         2.0.1     2021-08-05 [1] CRAN (R 4.1.0)
#>  reshape2       1.4.4     2020-04-09 [1] CRAN (R 4.1.0)
#>  rlang          1.0.4     2022-07-12 [1] CRAN (R 4.1.2)
#>  rmarkdown      2.14      2022-04-25 [1] CRAN (R 4.1.2)
#>  rncl           0.8.6     2022-03-18 [1] CRAN (R 4.1.2)
#>  RNeXML         2.4.7     2022-05-13 [1] CRAN (R 4.1.2)
#>  RSpectra       0.16-1    2022-04-24 [1] CRAN (R 4.1.2)
#>  rstudioapi     0.13      2020-11-12 [1] CRAN (R 4.1.0)
#>  scales         1.2.0     2022-04-13 [1] CRAN (R 4.1.2)
#>  sessioninfo    1.2.2     2021-12-06 [1] CRAN (R 4.1.0)
#>  stringi        1.7.8     2022-07-11 [1] CRAN (R 4.1.2)
#>  stringr        1.4.0     2019-02-10 [1] CRAN (R 4.1.0)
#>  styler         1.7.0     2022-03-13 [1] CRAN (R 4.1.2)
#>  tibble         3.1.8     2022-07-22 [1] CRAN (R 4.1.2)
#>  tidyr          1.2.0     2022-02-01 [1] CRAN (R 4.1.2)
#>  tidyselect     1.1.2     2022-02-21 [1] CRAN (R 4.1.2)
#>  tidytree       0.3.9     2022-03-04 [1] CRAN (R 4.1.2)
#>  treeio         1.21.0    2022-07-08 [1] Github (GuangchuangYu/treeio@9623fd6)
#>  utf8           1.2.2     2021-07-24 [1] CRAN (R 4.1.0)
#>  uuid           1.1-0     2022-04-19 [1] CRAN (R 4.1.2)
#>  vctrs          0.4.1     2022-04-13 [1] CRAN (R 4.1.2)
#>  withr          2.5.0     2022-03-03 [1] CRAN (R 4.1.2)
#>  xfun           0.31      2022-05-10 [1] CRAN (R 4.1.2)
#>  XML            3.99-0.10 2022-06-09 [1] CRAN (R 4.1.2)
#>  xml2           1.3.3     2021-11-30 [1] CRAN (R 4.1.0)
#>  yaml           2.3.5     2022-02-21 [1] CRAN (R 4.1.2)
#>  yulab.utils    0.0.5     2022-06-30 [1] CRAN (R 4.1.2)
#> 
#>  [1] /Library/Frameworks/R.framework/Versions/4.1/Resources/library
#> 
#> ──────────────────────────────────────────────────────────────────────────────

Release DAISIEprep 0.4.0

Prepare for release:

Submit to CRAN:

usethis::use_version('minor')
devtools::submit_cran()
Approve email

Wait for CRAN...

Accepted 🎉
~~Finish & publish blog post~~
~~Add link to blog post in pkgdown news menu~~
usethis::use_github_release()
~~usethis::use_dev_version(push = TRUE)~~ (at the moment not using development versioning)
~~Tweet~~

Add package documentation file

Would be useful to have documentation when running ?packagename for DAISIEprep.

Finish ASR vignette BioGeoBEARS example

`sse_states_to_endemicity()` and `endemicity_to_sse_states()` are non-symmetrical

I would expect sse_states_to_endemicity() and endemicity_to_sse_states() to reverse each others output. However, this is not the case.

endemicity_to_sse_states(endemicity_status = "nonendemic", sse_model = "geosse")

gives 0.

However

sse_states_to_endemicity(states = 0, sse_model = "geosse")

gives an error and says states should be 1, 2 or 3.

sse_states_to_endemicity(states = 3, sse_model = "geosse")

gives "0" (as a character) which is weird behaviour.

Perhaps all of this behaviour is intended. If so just let me know and we can document the functions to reflect that and I'll close the issue.

Draft ASR vignette introduction section

This section should explain the reader why it is relevant to perform ASR for DAISIE data preparation.
It should also describe an overview of the examples that are going to be covered through the vignette.

`add_asr_node_states` does not work with Mk model yet

From the Introduction section of the ASR vignette:

library(DAISIEprep)
set.seed(
  1,
  kind = "Mersenne-Twister",
  normal.kind = "Inversion",
  sample.kind = "Rejection"
)
phylo <- ape::rcoal(10)

phylo$tip.label <- c("bird_a", "bird_b", "bird_c", "bird_d", "bird_e", "bird_f",
                     "bird_g", "bird_h", "bird_i", "bird_j")

phylo <- phylobase::phylo4(phylo)

endemicity_status <- sample(
  x = c("not_present", "endemic", "nonendemic"),
  size = length(phylobase::tipLabels(phylo)),
  replace = TRUE,
  prob = c(0.6, 0.2, 0.2)
)

phylod <- phylobase::phylo4d(phylo, as.data.frame(endemicity_status))

# reconstruction using parsimony
phylod_parsimony <- add_asr_node_states(
  phylod = phylod,
  asr_method = "parsimony")

# reconstruction using Mk model
phylod_parsimony <- add_asr_node_states(
  phylod = phylod,
  asr_method = "Mk"
)
#> Error in ncol(asr$ancestral_likelihoods): object 'asr' not found

^{Created on 2022-07-28 by the reprex package (v2.0.1)}