Effect of character correlation

To measure the effect of the our inapplicable algorithm, compared to treating inapplicable tokens as missing data or an extra state, we ran a normal parsimony search on 22 matrices from 10 taxonomic classes (plants and animals) with variable number of characters, taxa and inapplicable data.

group taxa characters % NA characters % NA cells
Aguado2009 Annelids 76 107 40.2 5.0
Asher2005 Mammals 23 126 31.7 3.9
Baker2009 Angiosperms 205 105 35.2 8.3
Bouchenak2010 Angiosperms 818 417 81.5 21.1
Dikow2009 Insects 88 220 8.2 1.6
Eklund2004 Angiosperms 54 131 31.3 10.0
Geisler2001 Mammals 68 186 46.8 2.0
Giles2015 Fishes 78 236 74.2 17.3
Griswold1999 Spiders 43 137 43.8 9.6
Liljeblad2008 Insects 68 308 31.5 6.4
Loconte1991 Angiosperms 56 104 63.5 4.1
Longrich2010 Dinosaurs 19 119 26.1 5.7
OMeara2013 Mammals 63 317 44.5 5.4
Rougier2012 Mammals 58 317 2.8 0.4
Sharkey2011 Insects 105 392 100.0 20.1
Sundue2010 Pteridophytes 177 237 81.9 26.5
Vinther2008 Annelids 23 57 43.9 20.9
Wiens2010 Squamates 64 363 41.9 7.7
Wilson2003 Crustaceans 61 165 43.0 11.1
Wortley2006 Angiosperms 37 105 20.0 3.7
Zanol2012 Annelids 80 230 69.1 17.3
Zhu2013 Fishes 75 253 70.8 19.3

From these matrices, we ran a heuristic search in PAUP* treating all inapplicable data as missing, and saving the island of n trees (default PAUP* options). The length of of the trees on this island is called SMP the overall most parsimonious score. All trees on this island have the same length SMP.

The length of the trees from these island are then re-calculated by treating inapplicable as an extra state (or “new state”) or by using our inapplicable algorithm. This results in two new lists of n scores (for the n trees of each island) with SNS and SNA being the variable tree lengths for each topology under respectively the “new state” and inapplicable counting algorithms.

We then calculated how many trees were not the shortest in SNS and SNA (i.e. how many scores in these distributions were greater than the minimum) to measure how many trees were not discarded when treating inapplicable as missing data. Finally, we calculated how longer (or shorter!) were the SNS and SNA scores compared to the SMP to see how both “new state” and the inapplicable algorithm modify tree length in comparison to treating inapplicable data as missing data.

Before starting

Loading the functions

Getting the functions

source("Functions/read.tree.score.R")
source("Functions/standardise.score.R")
source("Functions/plot.scores.R")
source("Functions/test.R")

Loading the data

The most parsimonious islands treating inapplicable data as missing were obtained using PAUP* (v 4.0a151) with random sequences addition replicated 100 times with a maximum of 5 x 10^6 rearrangements per replicates (hsearch addseq=random nreps=100 rearrlimit=5000000 limitperrep=yes in PAUP). The length of the resulting trees within the islands were then measured using PAWM with inapplicable data treated as a new state or using our inapplicable algorithm.

We can load these results using the following:

## Setting the path
path <- "Data/Scores/"

## Getting the chain names
chains <- list.files(path, pattern = ".morphy.inapplicable.count")
chains <- unlist(strsplit(chains, split = ".morphy.inapplicable.count"))

## Remove Sutton2012 (no NA)
chains <- chains[-which(chains == "Sutton2012")]

## Loading the tree scores
scores <- sapply(chains, read.tree.score, path, simplify = FALSE)

## Loading the matrices
matrices <- lapply(chains, read.matrices, path = "Data/Matrices")

The scores are then standardise to reflect the proportional length of the trees compared to the shortest tree for both SNS and SNA.

## Normalising the score (dividing by the shortest score from each individual search)
scores_short <- mapply(standardise.score, scores, std = "short", SIMPLIFY = FALSE)

Getting the proportion of discarded trees under both algorithms treating inapplicable data as non missing

We can then measure from these standardised scores, how many trees would have been discarded from a default PAUP search (i.e. how values in SNS and SNA are higher than their minimum). Here, the proportion of discarded trees are the number of trees in the most parsimonious island (under each algorithm) that are one step or more longer than the shortest tree.

## Getting the proportion of discarded trees
proportions <- lapply(scores_short, get.proportion)
## Only selecting these proportions for inapplicable or extra state
inapp_proportions <- unlist(lapply(proportions, `[[`, 1))
extra_proportions <- unlist(lapply(proportions, `[[`, 3))
## Combining the data into a data frame
proportions_combined <- data.frame(inapp_proportions, extra_proportions)
colnames(proportions_combined) <- c("Inapplicable", "New state")

We can then visualise this as follows:

op <- par(bty = "n")
boxplot(proportions_combined, ylab = "Proportion of discarded (longer) trees")

par(op)
Min. 1st Qu. Median Mean 3rd Qu. Max.
Inapplicable 0 0.38 0.69 0.62 0.99 1
New state 0 0.50 0.76 0.70 0.99 1

We see here a clear effect of the counting algorithm used with nearly 100% of trees discarded (the results are rounded) when not treating inapplicable data as missing data. All these trees would have still be saved in a default analysis. Note that, although their quantiles overlap, our inapplicable algorithm produces, on average, slightly less trees that are longer than treating inapplicable data as a new state.

How long are the trees?

This results however, only indicate how many longer trees are generated under these two counting methods, not how long each of them are. We can compare each tree length SNS and SNA to their length SMP when treating inapplicable data as missing. This information can be added to the previous graph where in a two dimensional way: the vertical information is the same as above (proportion of discarded trees) and the horizontal gives information on, when trees are discarded, what are their minimum and maximum length compared to SMP.

## Normalising the score (based on the MP island)
scores_MP <- mapply(standardise.score, scores, std = "MP", SIMPLIFY = FALSE)
## Getting the proportion of discarded trees
proportions_MP <- lapply(scores_MP, get.proportion)
## Only selecting these proportions for inapplicable or extra state
inapp_proportions_MP <- unlist(lapply(proportions_MP, `[[`, 1))
extra_proportions_MP <- unlist(lapply(proportions_MP, `[[`, 3))
## Combining the data into a data frame
proportions_MP <- data.frame(inapp_proportions_MP, extra_proportions_MP)
colnames(proportions_MP) <- c("Inapplicable", "New state")
## Getting the relative minimum and maximum tree length from both distributions
max_scores <- lapply(scores_MP, lapply, max)
min_scores <- lapply(scores_MP, lapply, min)
## Only selecting these proportions for inapplicable or extra state
inapp_max <- unlist(lapply(max_scores, `[[`, 1))
extra_max <- unlist(lapply(max_scores, `[[`, 3))
inapp_min <- unlist(lapply(min_scores, `[[`, 1))
extra_min <- unlist(lapply(min_scores, `[[`, 3))
## Updating the data frame (using the inapp and extra proportion from before)
proportions_combined <- data.frame(inapp_proportions, inapp_min, inapp_max,
    extra_proportions, extra_min, extra_max)
op <- par(bty = "n")
sauronplot(proportions_combined, names = c("Inapplicable", "New state"),
    xlab = "Max/min proportional extra length to MP trees",
    ylab = "Proportion of discarded trees")

par(op)

Again, on the y axis is the same as above but represented as a line: the black dot is the median and the solid and dashed lines are the 50% and 95% confidence intervals respectively. The relative additional tree length to the overall most parsimonious tree is represented on the x axis. Each deviation represents the proportional extra length of the shortest and longest tree respectively on the left and right from the 0 line.

Min. 1st Qu. Median Mean 3rd Qu. Max.
Minimum extra length (NA) 0.00 0.01 0.04 0.10 0.10 0.87
Maximum extra length (NA) 0.00 0.01 0.05 0.12 0.11 0.90
Minimum extra length (NS) 0.01 0.18 0.25 0.56 0.60 3.10
Maximum extra length (NS) 0.01 0.18 0.26 0.58 0.65 3.12

Although both algorithms can generate longer trees (even for the minimum extra length), our inapplicable algorithm always generates trees less than 50% longer than the overall most parsimonious trees (apart from some trees were nearly every trees would have been discarded in the search).