---
title: "24_Nov_Oro_BAMM"
author: "Sebastian"
date: "11/24/2018"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```


Load required package and set wd
```{r}

require(ape)
require(BAMMtools)
library(coda)
setwd("~/Documents/pyphlawd_local/oro_bamm/24NOV18")


```


read in tree:

```{r}

tree <- read.tree("Fu_oro_multiCalib.summed.dated.newick.tre")

```


get data from each run, remove 10% as burn in, and check out convergence
```{r}

edata0.1 <- getEventData(tree, eventdata = "./24_nov_bamm/0.1/oro_bamm_24nov_r0.1_event_data.txt", burnin=0.1)
edata1.0 <- getEventData(tree, eventdata = "./24_nov_bamm/1.0/oro_bamm_24nov_r1.0_event_data.txt", burnin=0.1)
edata1.5 <- getEventData(tree, eventdata = "./24_nov_bamm/1.5/oro_bamm_24nov_r1.5_event_data.txt", burnin=0.1)
edata2.0 <- getEventData(tree, eventdata = "./24_nov_bamm/2.0/oro_bamm_24nov_r2.0_event_data.txt", burnin=0.1)
edata10.0 <- getEventData(tree, eventdata = "./24_nov_bamm/10.0/oro_bamm_24nov_r10.0_event_data.txt", burnin=0.1)
edata25.0 <- getEventData(tree, eventdata = "./24_nov_bamm/25.0/oro_bamm_24nov_r25.0_event_data.txt", burnin=0.1)

mcmcout0.1 <- read.csv("./24_nov_bamm/0.1/oro_bamm_24nov_r0.1_mcmc_out.txt", header=T)
mcmcout1.0 <- read.csv("./24_nov_bamm/1.0/oro_bamm_24nov_r1.0_mcmc_out.txt", header=T)
mcmcout1.5 <- read.csv("./24_nov_bamm/1.5/oro_bamm_24nov_r1.5_mcmc_out.txt", header=T)
mcmcout2.0 <- read.csv("./24_nov_bamm/2.0/oro_bamm_24nov_r2.0_mcmc_out.txt", header=T)
mcmcout10.0 <- read.csv("./24_nov_bamm/10.0/oro_bamm_24nov_r10.0_mcmc_out.txt", header=T)
mcmcout25.0 <- read.csv("./24_nov_bamm/25.0/oro_bamm_24nov_r25.0_mcmc_out.txt", header=T)


#mcmc
pdf("mcmc_likelihoods.pdf",width=10,height=10)
par(mfrow=c(2,3), mar=c(4.5, 2.0, 2.0, 1.0), xpd=TRUE)
plot(mcmcout0.1$logLik ~ mcmcout0.1$generation)
title(main="0.1")
plot(mcmcout1.0$logLik ~ mcmcout1.0$generation) #looks ok
title(main="1.0")
plot(mcmcout1.5$logLik ~ mcmcout1.0$generation) #looks ok
title(main="1.5")
plot(mcmcout1.0$logLik ~ mcmcout1.0$generation) #looks ok
title(main="2.0")
plot(mcmcout1.0$logLik ~ mcmcout1.0$generation) #looks ok
title(main="10.0")
plot(mcmcout1.0$logLik ~ mcmcout1.0$generation) #looks ok
title(main="25.0")
dev.off()

```


ESS with first 10% of samples
```{r}

#ESS
burnstart0.1 <- floor(0.1 * nrow(mcmcout0.1))
postburn0.1 <- mcmcout0.1[burnstart0.1:nrow(mcmcout0.1), ]
effectiveSize(postburn0.1$N_shifts)
effectiveSize(postburn0.1$logLik)

burnstart1.0 <- floor(0.1 * nrow(mcmcout1.0))
postburn1.0 <- mcmcout1.0[burnstart1.0:nrow(mcmcout1.0), ]
effectiveSize(postburn1.0$N_shifts)
effectiveSize(postburn1.0$logLik)

burnstart1.5 <- floor(0.1 * nrow(mcmcout1.5))
postburn1.5 <- mcmcout1.5[burnstart1.5:nrow(mcmcout1.5), ]
effectiveSize(postburn1.5$N_shifts)
effectiveSize(postburn1.5$logLik)

burnstart2.0 <- floor(0.1 * nrow(mcmcout2.0))
postburn2.0 <- mcmcout2.0[burnstart2.0:nrow(mcmcout2.0), ]
effectiveSize(postburn2.0$N_shifts)
effectiveSize(postburn2.0$logLik)

burnstart10.0 <- floor(0.1 * nrow(mcmcout10.0))
postburn10.0 <- mcmcout10.0[burnstart10.0:nrow(mcmcout10.0), ]
effectiveSize(postburn10.0$N_shifts)
effectiveSize(postburn10.0$logLik)

burnstart25.0 <- floor(0.1 * nrow(mcmcout25.0))
postburn25.0 <- mcmcout25.0[burnstart25.0:nrow(mcmcout25.0), ]
effectiveSize(postburn25.0$N_shifts)
effectiveSize(postburn25.0$logLik)


```


Prior vs. Posterior plots
```{r}

pdf("post_vs_prior.pdf",width=10,height=10)
par(mfrow=c(2,3), mar=c(3.0, 4.0, 1.5, 1.5), xpd=TRUE, cex=0.4)
plotPrior(mcmcout0.1, expectedNumberOfShifts=0.1, priorCol = 'grey', postCol = 'red')
title(main="0.1")
plotPrior(mcmcout1.0, expectedNumberOfShifts=1.0, priorCol = 'grey', postCol = 'red')
title(main="1.0")
plotPrior(mcmcout1.5, expectedNumberOfShifts=1.5, priorCol = 'grey', postCol = 'red')
title(main="1.5")
plotPrior(mcmcout2.0, expectedNumberOfShifts=2.0, priorCol = 'grey', postCol = 'red')
title(main="2.0")
plotPrior(mcmcout10.0, expectedNumberOfShifts=10.0, priorCol = 'grey', postCol = 'red')
title(main="10.0")
plotPrior(mcmcout25.0, expectedNumberOfShifts=25.0, priorCol = 'grey', postCol = 'red')
title(main="25.0")
dev.off()

```



distinct shift configurations:
```{r}

css0.1 <- credibleShiftSet(edata0.1, expectedNumberOfShifts=0.1, threshold=5, set.limit = 0.95)
css1.0 <- credibleShiftSet(edata1.0, expectedNumberOfShifts=1.0, threshold=5, set.limit = 0.95)
css1.5 <- credibleShiftSet(edata1.5, expectedNumberOfShifts=1.5, threshold=5, set.limit = 0.95)
css2.0 <- credibleShiftSet(edata2.0, expectedNumberOfShifts=2.0, threshold=5, set.limit = 0.95)
css10.0 <- credibleShiftSet(edata10.0, expectedNumberOfShifts=10.0, threshold=5, set.limit = 0.95)
css25.0 <- credibleShiftSet(edata25.0, expectedNumberOfShifts=25.0, threshold=5, set.limit = 0.95)


css0.1$number.distinct
css1.0$number.distinct
css1.5$number.distinct
css2.0$number.distinct
css10.0$number.distinct
css25.0$number.distinct

par(mfrow=c(2,3), mar=c(1, 0.5, 0.5, 0.5), xpd=TRUE)
plot.credibleshiftset(css0.1)
plot.credibleshiftset(css1.0)
plot.credibleshiftset(css1.5)
plot.credibleshiftset(css2.0)
plot.credibleshiftset(css10.0)
plot.credibleshiftset(css25.0)

```


Best shift configuration for each analysis
```{r}

best0.1 <- getBestShiftConfiguration(edata0.1, expectedNumberOfShifts=0.1)
best1.0 <- getBestShiftConfiguration(edata1.0, expectedNumberOfShifts=1.0)
best1.5 <- getBestShiftConfiguration(edata1.5, expectedNumberOfShifts=1.5)
best2.0 <- getBestShiftConfiguration(edata2.0, expectedNumberOfShifts=2.0)
best10.0 <- getBestShiftConfiguration(edata10.0, expectedNumberOfShifts=10.0)
best25.0 <- getBestShiftConfiguration(edata25.0, expectedNumberOfShifts=25.0)

pdf("best_shift_configs_withLegend2.pdf",width=20,height=25)
par(mfrow=c(2,3), mar=c(3.0, 4.0, 1.5, 1.5), xpd=TRUE, cex=0.8)
q<-plot.bammdata(best0.1, lwd = 1.5, legend = T, cex = 0.1)
title(main="0.1")
addBAMMshifts(best0.1, cex=2.0, col = 'red', par.reset = FALSE)
q<-plot.bammdata(best1.0, lwd = 1.5, legend = T, cex = 0.1)
title(main="1.0")
addBAMMshifts(best1.0, cex=2.0, col = 'red', par.reset = FALSE )
q<-plot.bammdata(best1.5, lwd = 1.5, legend = T, cex = 0.1)
title(main="1.5")
addBAMMshifts(best1.5, cex=2.0, col = 'red', par.reset = FALSE )
q<-plot.bammdata(best2.0, lwd = 1.5, legend = T, cex = 0.1)
title(main="2.0")
addBAMMshifts(best2.0, cex=2.0, col = 'red', par.reset = FALSE )
q<-plot.bammdata(best10.0, lwd = 1.5, legend = T, cex = 0.1)
title(main="10.0")
addBAMMshifts(best10.0, cex=2.0, col = 'red', par.reset = FALSE )
q<-plot.bammdata(best25.0, lwd = 1.5, legend = T, cex = 0.1)
title(main="25.0")
addBAMMshifts(best25.0, cex=2.0, col = 'red', par.reset = FALSE )
dev.off()


```


rate through time plot
```{r}

pdf("rate_through_time.pdf",width=10,height=10)
par(mfrow=c(1,1), mar=c(3.0, 4.0, 1.5, 1.5), xpd=TRUE, cex=0.4)
plotRateThroughTime(edata0.1, ratetype="speciation", xlim = c(0, 45), ylim = c(0, 0.6), avgCol = 'white', intervalCol = 'lightblue')
par(new=T)
plotRateThroughTime(edata1.0, ratetype="speciation", xlim = c(0, 45), ylim = c(0, 0.6), avgCol = 'white', intervalCol = 'lightblue')
par(new=T)
plotRateThroughTime(edata1.5, ratetype="speciation", xlim = c(0, 45), ylim = c(0, 0.6), avgCol = 'white', intervalCol = 'lightblue')
par(new=T)
plotRateThroughTime(edata2.0, ratetype="speciation", xlim = c(0, 45), ylim = c(0, 0.6), avgCol = 'white', intervalCol = 'lightblue')
par(new=T)
plotRateThroughTime(edata10.0, ratetype="speciation", xlim = c(0, 45), ylim = c(0, 0.6), avgCol = 'white', intervalCol = 'lightblue')
par(new=T)
plotRateThroughTime(edata25.0, ratetype="speciation", xlim = c(0, 45), ylim = c(0, 0.6), avgCol = 'white', intervalCol = 'lightblue')
par(new=T)
plotRateThroughTime(edata0.1, ratetype="speciation", xlim = c(0, 45), ylim = c(0, 0.6), avgCol = 'red', intervalCol = NULL)
par(new=T)
plotRateThroughTime(edata1.0, ratetype="speciation", xlim = c(0, 45), ylim = c(0, 0.6), avgCol = 'orange', intervalCol = NULL)
par(new=T)
plotRateThroughTime(edata1.5, ratetype="speciation", xlim = c(0, 45), ylim = c(0, 0.6), avgCol = 'yellow', intervalCol = NULL)
par(new=T)
plotRateThroughTime(edata2.0, ratetype="speciation", xlim = c(0, 45), ylim = c(0, 0.6), avgCol = 'green', intervalCol = NULL)
par(new=T)
plotRateThroughTime(edata10.0, ratetype="speciation", xlim = c(0, 45), ylim = c(0, 0.6), avgCol = 'blue', intervalCol = NULL)
par(new=T)
plotRateThroughTime(edata25.0, ratetype="speciation", xlim = c(0, 45), ylim = c(0, 0.6), avgCol = 'violet', intervalCol = NULL)
dev.off()

```


plots of each best shift configuration with tip labels
```{r}



pdf("best_shift_config_0.1.pdf",width=10,height=40)
q<-plot.bammdata(best0.1, lwd = 1.0, labels = TRUE, cex = 0.25)
title(main="0.1")
addBAMMshifts(best0.1, cex=1.5, col = 'pink', par.reset = FALSE )
dev.off()


pdf("best_shift_config_1.0.pdf",width=10,height=40)
q<-plot.bammdata(best1.0, lwd = 1.0, labels = TRUE, cex = 0.25)
title(main="1.0")
addBAMMshifts(best1.0, cex=1.5, col = 'pink', par.reset = FALSE )


pdf("best_shift_config_1.5.pdf",width=10,height=40)
q<-plot.bammdata(best1.5, lwd = 1.0, labels = TRUE, cex = 0.25)
title(main="1.5")
addBAMMshifts(best1.5, cex=1.5, col = 'pink', par.reset = FALSE )


pdf("best_shift_config_2.0.pdf",width=10,height=40)
q<-plot.bammdata(best2.0, lwd = 1.0, labels = TRUE, cex = 0.25)
title(main="2.0")
addBAMMshifts(best2.0, cex=1.5, col = 'pink', par.reset = FALSE )


pdf("best_shift_config_10.0.pdf",width=10,height=40)
q<-plot.bammdata(best10.0, lwd = 1.0, labels = TRUE, cex = 0.25)
title(main="10.0")
addBAMMshifts(best10.0, cex=1.5, col = 'pink', par.reset = FALSE )


pdf("best_shift_config_25.0.pdf",width=10,height=40)
q<-plot.bammdata(best25.0, lwd = 1.0, labels = TRUE, cex = 0.25)
title(main="25.0")
addBAMMshifts(best25.0, cex=1.5, col = 'pink', par.reset = FALSE )
dev.off()




```


```{r}

best0.1$numberEvents
best1.0$numberEvents
best1.5$numberEvents
best2.0$numberEvents
best10.0$numberEvents
best25.0$numberEvents



```

```{r}

pdf("best_shift_config_1.0_forPresentation.pdf",width=30,height=80)
q<-plot.bammdata(best1.0, lwd = 15.0, labels = TRUE, cex = 0.25)
title(main="1.0")
addBAMMshifts(best1.0, cex=10, col = 'pink', par.reset = FALSE )






```

