---
title: "TreeSearch results"
author: "Martin R. Smith"
---

# Corrected parsimony {#treesearch}

The phylogenetic dataset contains a considerable proportion of inapplicable codings (`r sum(my_chars == '-')`, = `r round(100 * sum(my_chars == '-') / sum(my_chars != '?'), 1)`% of `r sum(my_chars != '?')` non-ambiguous tokens;  `r round(sum(my_chars == '-') / length(my_chars), 3) * 100`% of `r length(my_chars)` total cells), which are known to introduce error and bias to phylogenetic reconstruction when the Fitch algorithm is employed [@Maddison1993;@Brazeau2018].  As such, we used the R package _TreeSearch_ v0.1.2 [@Smith2018TreeSearch] to conduct phylogenetic tree search with a tree-scoring algorithm that avoids logically impossible character transformations when handling inapplicable data [@Brazeau2018], implemented in the _MorphyLib_ C library [@Brazeau2017Morphylib].

## Search parameters

Heuristic searches were conducted using the parsimony ratchet [@Nixon1999] under equal and implied weights [@Goloboff1997].  The consensus tree presented in the main manuscript represents a strict consensus of all trees that are most parsimonious under one or more of the concavity constants (_k_) `r paste(kValues[-length(kValues)], collapse=', ')` and `r kValues[length(kValues)]`, an approach that has been shown to produce higher accuracy (i.e. more nodes and quartets resolved correctly) than equal weights at any given level of precision [@Smith2019].

## Analysis 

The R commands used to conduct the analysis are reproduced below.
The results can most readily be replicated using the [R markdown files](https://github.com/ms609/hyoliths/) used to generate these pages:
in Rstudio, run `r GitLink("index.Rmd", "index.Rmd")`, then run each block in 
`r GitLink("13_TreeSearch.Rmd", "TreeSearch.Rmd")`.  The complete analysis will take several hours.

### Initialize and load data

```{r treesearch-load-morphoBank, echo=TRUE, eval=FALSE}
# Load data from locally downloaded copy of MorphoBank matrix
my_data <- ReadAsPhyDat(nexusFile)
my_data[ignoredTaxa] <- NULL
iw_data <- PrepareDataIW(my_data)
```

### Generate starting tree

Start by quickly rearranging a neighbour-joining tree, rooted on the outgroup.

```{R treesearch-starting-tree, echo=TRUE}
nj.tree <- NJTree(my_data)
rooted.tree <- EnforceOutgroup(nj.tree, outgroup)
start.tree <- TreeSearch(tree=rooted.tree, dataset=my_data, maxIter=3000,
                         EdgeSwapper=RootedNNISwap, verbosity=0)
```

### Implied weights analysis

The position of the root does not affect tree score, so we keep it fixed (using `RootedXXXSwap` functions) to avoid unnecessary swaps.

```{r treesearch-implied-weights-analysis, echo=TRUE, eval=FALSE}
for (k in kValues) {
  iw.tree <- IWRatchet(start.tree, iw_data, concavity=k,
                       ratchHits = 20, ratchIter=4000, searchHits=56,
                       swappers=list(RootedTBRSwap, RootedSPRSwap, RootedNNISwap),
                       verbosity=0L)
  score <- IWScore(iw.tree, iw_data, concavity=k)
  # Write a single best tree
  write.nexus(iw.tree,
              file=paste0("TreeSearch/hk_iw_k", k, "_", 
                          signif(score, 5), ".nex", collapse=''))

  iw.consensus <- IWRatchetConsensus(iw.tree, iw_data, concavity=k,
                  swappers=list(RootedTBRSwap, RootedNNISwap),
                  searchHits=55, searchIter=4000, nSearch=250, verbosity=0L)
  write.nexus(iw.consensus,
              file=paste0("TreeSearch/hk_iw_k", k, "_", 
                          signif(IWScore(iw.consensus[[1]], iw_data, concavity=k), 5),
                          ".all.nex", collapse=''))
}
```


```{R treesearch-deeper-iw, echo=FALSE, eval=FALSE}
# This block provides a crude check that each run has found an optimum;
# if not it continues searching, and you may consider increasing the 
# ratchHits, ratchIter or searchHits parameters above.
for (k in kValues) {
  source('loadTrees.R')
  
  k.scores <- unlist(lapply(iw.exist, function (x) IWScore(x[[1]], iw_data, k)))
  ew.scores <- vapply(ew.trees, IWScore, double(1), iw_data, k)
  this.k <- kValues[iw.treesLoaded] == k
  if (!any(this.k)) next
  my.score <- k.scores[this.k]
  ew.best <- min(ew.scores)
  if (ew.best < min(k.scores)) {
    start.tree <- ew.trees[which.min(ew.scores)]
  } else if (k.scores[this.k] > min(k.scores)) {
    start.tree <- iw.trees[[which.min(k.scores)]][[1]]
  } else {
    cat("\n\n * The best tree at k =", k, "is better than any other optimal tree.")
    next
  }
  
  attr(start.tree, 'score') <- NULL
  cat("\n\n * Did not find global optimum for k =", k, "; continuing search.\n   ")
  iw.tree <- IWRatchet(start.tree, iw_data, concavity=k,
                       ratchHits = 30, ratchIter=4000, searchHits=56,
                       swappers=list(RootedTBRSwap, RootedSPRSwap, RootedNNISwap),
                       verbosity=1L)
  score <- IWScore(iw.tree, iw_data, concavity=k)
  # Write a single best tree
  write.nexus(iw.tree,
              file=paste0("TreeSearch/hk_iw_k", k, "_", 
                          signif(score, 5), ".nex", collapse=''))

  iw.consensus <- IWRatchetConsensus(iw.tree, iw_data, concavity=k,
                  swappers=list(RootedTBRSwap, RootedNNISwap),
                  searchHits=55, nSearch=250, searchIter=4000, verbosity=1L)
  write.nexus(iw.consensus,
              file=paste0("TreeSearch/hk_iw_k", k, "_", 
                          signif(IWScore(iw.consensus[[1]], iw_data, concavity=k), 5),
                          ".all.nex", collapse=''))
}
```

### Equal weights analysis

```{r treesearch-equal-weights-analysis, echo=TRUE, eval=FALSE}
ew.tree <- Ratchet(start.tree, my_data, verbosity=0L,
                   ratchHits = 20, ratchIter=4000, searchHits=55,
                   swappers=list(RootedTBRSwap, RootedSPRSwap, RootedNNISwap))
ew.consensus <- RatchetConsensus(ew.tree, my_data, nSearch=250, searchHits = 85,
                                 swappers=list(RootedTBRSwap, RootedNNISwap),
                                 verbosity=0L)
write.nexus(ew.consensus, file=paste0(collapse='', "TreeSearch/hk_ew_",
                                      Fitch(ew.tree, my_data), ".nex"))
```

## Results

Optimal trees can be downloaded in Nexus format from `r GitLink('TreeSearch', raw=FALSE)`.

```{R treesearch-load-trees, echo=FALSE}
source('loadTrees.R')
tipIndex <- sort(allTrees[[1]]$tip.label)
```

(ref:ts-all-cons) Consensus of all parsimony trees, under equal and implied weights.  

(ref:node-support) Node labels denote the frequency of each clade in most parsimonious trees under all analytical conditions.  

```{r treesearch-maj-consensus, fig.cap="(ref:ts-all-cons) (ref:node-support)", fig.width=6, fig.height=5.1, echo=FALSE}
majCon <- RootTree(consensus(allTrees, p=0.5), rootingTips)
supporters <- SplitFrequency(majCon, allTrees)

ColPlot(majCon)
nodelabels(paste0("\n\n", signif(supporters / length(allTrees), 2)),
           col=SupportColour(round(supporters / length(allTrees), 2)),
           adj=0, pos=2, frame='none', cex=0.75)

```

(ref:ts-trim-cons) Consensus of same trees, with taxa pruned before constructing consensus to give context to clade support.  

```{r treesearch-maj-consensus-pruned, fig.cap="(ref:ts-trim-cons) (ref:node-support)", fig.width=6, fig.height=5.1, echo=FALSE}
#omit <- c("Micrina", 'Mickwitzia_muralensis', 'Heliomedusa_orienta')
omit <- c("Cotyledion_tylodes", "Namacalathus", "Phoronis")
allPruned <- lapply(allTrees, drop.tip, omit)
majCon <- RootTree(consensus(allPruned, p=0.5), rootingTips)

####
supporters <- SplitFrequency(majCon, allPruned)
####
ColPlot(majCon)
nodelabels(paste0("\n", signif(supporters / length(allTrees), 2)),
           col=SupportColour(round(supporters / length(allTrees), 2)),
           adj=0, pos=2, frame='none', cex=0.75)
ColMissing(omit)

```

(ref:ts-iw-indiv) Strict consensus trees of implied weights analyses

\clearpage 

```{r treesearch-iw-indiv-1, fig.width=7.1, fig.height=5.1, fig.cap="(ref:ts-iw-indiv) (ref:first-panels)", echo=FALSE}
# Plot consensus results
par(mfrow=c(1, 2), mar=rep(0.2, 4))
PlotPanel(iw.trees, 1)
PlotPanel(iw.trees, 2)
```

\clearpage 

```{R treesearch-iw-indiv-2, fig.width=7.1, fig.height=5.1, fig.cap="(ref:ts-iw-indiv) (ref:second-panels)", echo=FALSE}
par(mfrow=c(1, 2), mar=rep(0.2, 4))
PlotPanel(iw.trees, 3)
PlotPanel(iw.trees, 4)
```

\clearpage 

```{R treesearch-iw-indiv-3, fig.width=7.1, fig.height=5.1, fig.cap="(ref:ts-iw-indiv) (ref:third-panels)", echo=FALSE}
par(mfrow=c(1, 2), mar=rep(0.2, 4))
PlotPanel(iw.trees, 5)
PlotPanel(iw.trees, 6)
```

\clearpage

```{r treesearch-equal-weights-results, fig.cap="Strict consensus of most parsimonious trees under equally weighted parsimony", fig.width=6, fig.height=7.1, echo=FALSE}
ew.trees <- lapply(ew.trees, drop.tip, 'Siphogonuchites_multa')
ColPlot(RootTree(consensus(ew.trees), rootingTips))
```


```{R treesearch-equal-weights-pruned-consensus, fig.cap="Strict consensus of equal weights results, taxa excluded", fig.width=6, fig.height=5.1, echo=FALSE, eval=FALSE}
omit <- c("Cotyledion_tylodes", "Micrina")
ColPlot(RootTree(ConsensusWithout(ew.trees, omit), rootingTips))
ColMissing(omit)
```

\clearpage
