wed_DEC_bsm.Rmd

---
output: html_document
---

```{r global_options, include=FALSE}
knitr::opts_chunk$set(fig.width=12, fig.height=8, eval = FALSE,
                      echo=TRUE, warning=FALSE, message=FALSE,
                      tidy = TRUE, collapse = TRUE,
                      results = 'hold')
```

In this exercise, you will use the results of the ancestral state recosntruction, to do a biogeographic stochastic mapping, and count the number of dispersals throught time.

# Library setup
```{r}
library(tidyverse)
library(ape)
library(geiger)
library(optimx)         # You need to have some version of optimx available
library(FD)       # for FD::maxent() (make sure this is up-to-date)
library(snow)     # (if you want to use multicore functionality; some systems/R versions prefer library(parallel), try either)
library(parallel)
library(devtools)
library(rexpokit)
library(cladoRcpp)
library(BioGeoBEARS)
library(stringr)
library(RColorBrewer)
library(colorspace)
library(jpeg)
library(viridis)
```



# Loading the results of the DEC ancestral area reconstruction

```{r}
load("example_data/bombacoideae_DEC_results.Rdata")
model_name = "DEC"
results_object = res = resDEC
scriptdir = np(system.file("extdata/a_scripts", package="BioGeoBEARS"))

clado_events_tables = NULL
ana_events_tables = NULL
lnum = 0

BSM_inputs_fn = "BSM_inputs_file.Rdata"
runInputsSlow = TRUE
if (runInputsSlow)
{
  stochastic_mapping_inputs_list = get_inputs_for_stochastic_mapping(res=res)
  save(stochastic_mapping_inputs_list, file=BSM_inputs_fn)
} else {
  # Loads to "stochastic_mapping_inputs_list"
  load(BSM_inputs_fn)
} # END if (runInputsSlow)
```

# Run the biogeographic stochastic mapping

```{r}
set.seed(seed=as.numeric(Sys.time()))

runBSMslow = TRUE

BSM_output = runBSM(res,
                    stochastic_mapping_inputs_list = stochastic_mapping_inputs_list, 
                    maxnum_maps_to_try=2000, 
                    nummaps_goal=10, # INcrease the number of replicates if possible
                    maxtries_per_branch=40000,
                    save_after_every_try=TRUE, 
                    savedir=getwd(), 
                    seedval=12345,
                    wait_before_save=0.01)

RES_clado_events_tables = BSM_output$RES_clado_events_tables
RES_ana_events_tables = BSM_output$RES_ana_events_tables
```

# Save the results to disk
```{r}
# Extract BSM output
clado_events_tables = BSM_output$RES_clado_events_tables
ana_events_tables = BSM_output$RES_ana_events_tables
save(clado_events_tables, file = "example_data/bombacoideae_BSM_clado_events_table_mcc_detbiom.rda")
save(ana_events_tables, file ="example_data/bombacoideae_BSMana_events_table_mcc_detbiom.rda")
```

## Calculate the number of dispersal through time
The following function calculates the number of estimated dispersals through time relative to the total branch length across time bins. 
The function requires the `BSM_output` and `res` objects created in the previous steps. 
Additionally you can adjust the max age (based on the root age of your tree) and the bin size. 
```{r}
count_dispersal_events <- function(BSM_output, res,
                                  max_age = 25,
                                  bin_size = 1,
                                  FUN = mean){
      resolution <- bin_size
      tr <- read.tree(res$inputs$trfn)
      
      clado_events_tables = BSM_output$RES_clado_events_tables
      ana_events_tables = BSM_output$RES_ana_events_tables
      
      area_labels  = sort(unique(ana_events_tables[[length(ana_events_tables)]]$dispersal_to))

      # Get the number of shifts through time
      breaks <- seq(max_age, 0, by = -resolution)
      labels <- rev((breaks+resolution/2)[-1])

      angen <- lapply(ana_events_tables, function(k) {
          out <- k %>% dplyr::select(abs_event_time, event_type, event_txt, dispersal_to) %>% mutate(bin = cut(abs_event_time, breaks = breaks, labels = labels)) %>% mutate(from = str_split_fixed(event_txt, n = 2, pattern = "->")[, 1]) %>% mutate(sdb_efb = ifelse(!grepl("W", from) & grepl("W", dispersal_to), 
              1, 0)) %>% mutate(efb_sdb = ifelse(from == "W" & (grepl("D", dispersal_to) | grepl("S", dispersal_to)), 1, 0))
          return(out)
      })


      # count events
      disps <- lapply(angen, function(k) {
          out <- k %>% filter(event_type == "d") %>% group_by(bin, dispersal_to) %>% summarize(abs_count = n()) %>% ungroup() %>% mutate(bin = parse_number(as.character(bin)))
    
          fram <- expand_grid(bin = labels, dispersal_to = area_labels)
          out <- left_join(fram, out, by = c("dispersal_to", "bin")) %>% replace_na(list(abs_count = 0))
      })

      disps <- bind_rows(disps, .id = "SM") %>% group_by(bin, dispersal_to) %>% summarize(lower = quantile(abs_count, probs = 0.025), avg = FUN(abs_count), upper = quantile(abs_count, 
          probs = 0.975))
      
      
      disps <- disps[which(disps[,5] > 0),]
      
      # get total branch length
      milgaps_mod <- matrix(ncol = 2, nrow = length(labels), 0)
      colnames(milgaps_mod) <- c("time", "tot_brl")
      milgaps_mod[, 1] <- labels
      rev_breaks = rev(breaks)
      ## Sum times per bin
      v <- branching.times(tr)
      for (i in 2:length(rev_breaks)) {
              s <- rev_breaks[i]  # time start
              e <- rev_breaks[i-1]  # time end
              l1 <- length(v[v > s]) * resolution + min(length(v[v > s]), 1) * resolution
              l2 <- sum(v[v > e & v < s] - e)
              if (l1 == 0) 
                  {
                      l2 = l2 * 2
                  }  # root
              milgaps_mod[i-1, 2] <- l1 + l2
      }
      brls <- data.frame(milgaps_mod)
      # 
      rescaled_disps <- left_join(disps, brls, by = c(bin = "time")) %>% mutate(lower = lower/tot_brl) %>% mutate(avg = avg/tot_brl) %>% mutate(upper = upper/tot_brl)
      rescaled_disps[,1] <- -rescaled_disps[,1]
      return(rescaled_disps)
    
}
```

After loading the function you can call it to get a table of relative dispersal counts and plot them through time
```{r}
relative_dispersals <- count_dispersal_events(BSM_output, res, max_age = 25, bin_size = 1)

# plot results
p <- ggplot() + geom_point(data = relative_dispersals, aes(x = bin, y = avg, color = dispersal_to)) + 
                geom_line(data = relative_dispersals, aes(x = bin, y = avg, color = dispersal_to)) + theme_bw()
p + labs(x = "Time (Ma)", y = "Relative number of dispersals")
```