---
title: "treestructure applied to structured coalescent simulation"
author: "Erik Volz"
date: "`r Sys.Date()`"
output:
  rmarkdown::html_vignette:
  #rmarkdown::html_vignette
  #bookdown::pdf_book:
    toc: TRUE
pkgdown:
  as_is: true
fontsize: 12pt
vignette: >
  %\VignetteIndexEntry{treestructure applied to structured coalescent simulation}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 11,
  warning = FALSE,
  message = FALSE
)
```


# Structured coalescent simulation

This example shows the function `trestruct` applied to a simulated structured
coalescent tree that includes samples from a large constant size population and
samples from three small "outbreaks" which are growing exponentially.
These simulations were generated with the [phydynR package](https://emvolz-phylodynamics.github.io/phydynR/).

```{r message=FALSE}
library(ape)
library(treestructure)
library(ggtree)
```

Load the tree:

```{r}
( tree <- ape::read.tree( system.file('sim.nwk', package = 'treestructure') ) )
```


Note that the tip labels corresponds to the deme of each sample.
'1' is the constant size reservoir, and '0' is the exponentially growing deme.


This will run the `treestructure` algorithm under default settings:
```{r message=FALSE}
s <- trestruct( tree )
```

The `treestructure` algorithm searches from root to tips of the phylogenetic tree.
The message printed on screen, `Finding splits under nodes`, reports which nodes
the algorithm is currently testing.

By default `trestruct` calibrates the split threshold to a target false discovery
rate (`fdr = 0.2`), so that the expected proportion of spuriously designated clusters
is controlled across the whole tree. The printed output reports this target and a
global test for whether the tree has any structure at all:

```{r}
print(s)
```

## Plotting results

The default plotting behavior uses the `ggtree` package if available.
```{r message=FALSE}
plot(s)  + ggtree::geom_tiplab()
```

If not, or if desired, `ape` plots are available
```{r}
plot( s, use_ggtree = FALSE )
```


For subsequent analysis, you may want to turn the `treestructure` result into a
dataframe:
```{r}
structureData <- as.data.frame( s )
head( structureData )
```

Each cluster and partition assignment is stored as a factor. You could use `split`
to get a data frame for each partition.
Suppose we want a tree corresponding to partition 1:
```{r}
with ( structureData,
       ape::keep.tip(s$tree, taxon[ partition==1 ] )
       ) -> partition1
partition1
plot(partition1)
```


# Parameter choice and number of clusters

Two things have the largest influence on the number of clusters.

1. `minCladeSize` controls the smallest allowed cluster in terms of the number of
tips; the default is 10. With a smaller value, smaller clusters may be detected, but
computation time increases. For example, structure that exists in clades smaller than
`minCladeSize` will not be detected regardless of the threshold, so lowering it can
reveal finer structure:

```{r message=FALSE}
trestruct( tree, minCladeSize = 5 )
```

2. The split threshold — how readily a clade is designated a new cluster. By default
this is set by a false discovery rate (`fdr`); alternatively it can be set by a
subjective significance `level`, or that `level` chosen automatically with the CH
index.

### False discovery rate (default)

By default the threshold is calibrated to a target false discovery rate, so that
`fdr` has an interpretable, whole-tree error-rate meaning:

```{r message=FALSE}
trestruct( tree, fdr = 0.05 )
```

The [false discovery rate](false_discovery_rate.html) vignette explains what the rate
controls, the global test for any structure, and the effect of serial (heterochronous)
sampling.

### Subjective significance level

Supplying a `level` (without `fdr`) uses it instead as a fixed per-test significance.
This is a subjective clustering threshold rather than an error rate: raising it detects
more clusters but also increases the false positive rate.

```{r message=FALSE}
trestruct( tree, level = 0.05 )
```

The best value of `level` depends on your application. One way to choose it is to use
additional data associated with each sample, selecting the `level` that gives clusters
explaining the most variance in a variable of interest (e.g. using the cluster as a
factor in an ANOVA).

### CH index

Alternatively, in the absence of any additional data, the `treestructure` package
supports using the [CH index](https://en.wikipedia.org/wiki/Calinski–Harabasz_index)
to compare different `level`s. This statistic is based on the ratio of the
between-cluster and within-cluster variance of the time of each node
(distance from the root) and returns the `level` such that this ratio is maximized.
If you wish to use the CH index, pass `level = NULL` to `trestruct`, and
read documentation for the `levellb`, `levelub`, and `res` parameters.