## TreePar {#treepar}

Let us walk here through a very simple TreePar analysis [@Stadler2011].
For a real publication worthy analysis, you may want to do a more thorough analysis by allowing more possible shifts and a finer grid for the search space.

First, load the TreePar package.
```{r load-TreePar}
  library(TreePar)
```
  
  TreePar needs the ages of the phylogeny to perform it's likelihood computation.
We use the TreePar builtin functionality for this.
```{r extract-ages, cache=TRUE}
x     <- sort(getx(tree),decreasing=TRUE)
max_t <- max(x)
```
For our toy analysis, we setup some simpler than normal search values.
We specify that there are at most three rate shifts, and that there are ten possible shift points.
For your own analysis, you should most likely increase these values.
We also specify an interval offset of 0.05, which means that in the most recent 5% of the tree age there was no rate shift event.
This is simply for improving the parameter estimation.
```{r TreePar-settings}
MAX.NUM.SHIFTS  = 4
NUM.INTERVALS   = 20
INTERVAL.OFFSET = 0.05
```
Now we run the full TreePar maximum likelihood search.
```{r TreePar-ML, echo=TRUE, cache=TRUE, results='hide'}
res   <- bd.shifts.optim(x,
                         sampling=rep(1,MAX.NUM.SHIFTS+1),
                         grid=max_t/NUM.INTERVALS,
                         start=INTERVAL.OFFSET*max_t,
                         end=max_t,
                         survival=1,
                         posdiv=TRUE)[[2]]
```
Now we have the maximum likelihood estimates for the different number of shifts.
We want to test if $i$ shifts explain the tree significantly better than $i-1$ shifts.
```{r TreePar-MLE}
best <- 0
for (i in 1:(length(res)-1)) {
    test <- pchisq(2*(res[[i]][1]-res[[i+1]][1]),3)
    #if test>0.95 then i shifts is significantly
    # better than i-1 shifts at a 5% error
    if ( test > 0.95 ) {
        best <- i
    } else {
        break
    }
}
```

```{r TreePar-rate-functions}
best_result <- res[[best+1]]
speciation <- c()
extinction <- c()
for (i in 1:(best+1)) {
    speciation[i] <- best_result[1+i+(best+1)] /
                     (1-best_result[1+i])
    extinction[i] <- best_result[1+i+(best+1)] /
                     (1/best_result[1+i]-1)
}

if ( best > 0 ) {
    shift_times <- best_result[ (length(best_result)-best+1) :
                                length(best_result)]
    tmp_lambda <- function(x) {
       index <- findInterval(x, shift_times);
       return ( speciation[index+1] )
    }
    tmp_mu     <- function(x) {
       index <- findInterval(x, shift_times);
       return ( extinction[index+1] )
    }
    lambda <- Vectorize(tmp_lambda)
    mu     <- Vectorize(tmp_mu)
} else {
    lambda <- function(x) { return ( rep(speciation[1],length(x)) ) }
    mu     <- function(x) { return ( rep(extinction[1],length(x)) ) }
}
```

Let us make a data frame representation of the rate functions
```{r treepar-rates-dataframe}
times <- seq(0, max_t, length.out = 100)

primates_ebd_treepar <- tibble(
  "time" = times,
  "lambda" = lambda(times),
  "mu" = mu(times),
)
head(primates_ebd_treepar)
```

Now let's have a look at the estimated rate functions
```{r plot-treepar-rates, fig.width=5, fig.height = 5, fig.cap = "The estimated posterior median results from an episodic birth-death model fitted in RevBayes."}

```

The data frame `primates_ebd_treepar` is what we get when we call the `data("primates_ebd_treepar")` function in `ACDC`. We use the RevBayes estimated rates in our vignette examples, but it is easy to substitute those for the two-epoch TreePar results to construct the congruence class.

