Analysis.qmd

---
title: "Analysis"
format: 
  html:
    embed-resources: true
    self-contained-math: true
editor: visual
editor_options: 
  chunk_output_type: console
---

```{r, echo=FALSE}
cacheData = F
```

# Read and process data

```{r}
library(RColorBrewer)
library(dplyr)
library(tidyr)
library(ranger)
library(stringr)
library(heatmaply)

# Read in simulated tree data
simuData <- read.csv("./experiments/uniform_sampling_experiment/process_parameter_values_and_tree_metrics_sign_corrected.csv", stringsAsFactors = T)
simuDataKey <- read.csv("./experiments/uniform_sampling_experiment/simulation_parameters_key.csv")

head(simuDataKey)

statisticsIndices1 = c(25, 27:43) # all 18 tree metrics
statisticsIndices2 = c(28, 29, 33:43) # # Only the tree metrics not strongly correlated with richness, --excluding log10S, PD, Colless, Sackin, Yule.PDA.ratio
statisticsIndices3 = c(25, 27:50) # all 18 tree metrics plus first 7 PCA

# Indices for the 5 simulation processes (env, dis, mut, nic, com)
predictorsIndices = c(14,16,18,20,22)
predictors = colnames(simuData)[c(14,16,18,20,22)]
numModelsPerPredictor = c(5, 7, 4, 7, 6)

# 8 models (accounting for abbreviation changes post-experiment)
modelAbbrevs = data.frame(original = c('ca', 'fh', 'gen', 'hs', 'la', 'pontarp', 've', 'xe'), new = c('ca', 'fh', 'gen', 'hs', 'am', 'pw', 've', 'xe'))

models = modelAbbrevs$original


```

# Clustering / redundancy analysis

TODO - add correlation analysis between metrics that was previously in the method section

```{r}
# Correlation among indices
corrs = cor(simuData[, statisticsIndices1], use = "pairwise.complete.obs", )

metricNames = colnames(simuData)[statisticsIndices1]
metricNames[metricNames == "mean.Iprime"] = "mean I'"

colnames(corrs) = metricNames
rownames(corrs) = metricNames

par(mar = c(8, 8, 10, 8))
heatmaply(corrs, colors = cool_warm)

```

# Consistency analysis

## Calculating correlation coefficients

Calculating individual as well as model-averaged correlation coefficients between all metrics and processes

```{r, cache=cacheData}
statisticsIndices = statisticsIndices1
statistics = colnames(simuData)[statisticsIndices]
nStats = length(statisticsIndices)

signMat = matrix(nrow = length(statistics), 
                 ncol = length(predictors),
                 dimnames = list(statistics, predictors))

R2Mat = disAg = disAg2 = signMat

R2MatFull = signMatFull = array(
    dim = c(length(statistics), length(predictors), length(models)),
    dimnames = list(statistics, predictors, models)) 


for(i in 1:length(statisticsIndices)){
  for(j in 1:length(predictorsIndices)){
    R2 = sig = rep(NA, length(models))
    for(k in 1:length(models)){
      tmp = simuData[simuData$model == models[k],]
      x = scale(tmp[,statisticsIndices[i]])
      y = tmp[,predictorsIndices[j]]
      if(! all(is.na(y))){
      
      corRes = cor.test(x,y, method = "kendall")
        
      # fit <- summary(lm(y ~ x))
      # R2[k] = fit$r.squared
      # sig[k] = ifelse(fit$coefficients[2,4] <0.05, sign(fit$coefficients[2,1]), 0)
      R2[k] = corRes$estimate
      sig[k] = corRes$p.value
      } else{
        R2[k] = NA
        sig[k] = NA
      }
    }
    R2MatFull[i,j,] = R2
    signMatFull[i,j,] = sig
    R2Mat[i,j] = mean(R2, na.rm = T)
    signMat[i,j] = sum(sig, na.rm = T) / sum(!is.na(sig))
    disAg[i,j] = abs(sum(sig == 1, na.rm = T) - sum(sig == -1, na.rm = T)) / abs(sum(sig %in% c(1,-1), na.rm = T))
    disAg2[i,j] = max(sum(sig == 1, na.rm = T), sum(sig == -1, na.rm = T)) / abs(sum(sig %in% c(1,-1), na.rm = T))
  }
}
```

## Full plot

Definition of image plot

```{r}
image.real3d <- function(mat, mat2, xCol = c("darkblue","blue", "lightblue", "gray94", "pink","red", "darkred"), 
                       range = c(-1,1), x.labels = dimnames(mat)[[1]],      y.labels = dimnames(mat)[[2]]) { 
  
  newMat = matrix(nrow = dim(mat)[1], ncol = 50)
  for(i in 1:dim(mat)[1]){
    for(k in 1:dim(mat)[2])
      newMat[i,(1+(k-1)*10):((k-1)*10 + dim(mat)[3])] = mat[i,k,]
  }
  
  newMat <- t(newMat)
  xpos = rep(0:(dim(mat)[2]-1), each = dim(mat)[3])*10 + 1:dim(mat)[3] 
  ypos = 1:dim(mat)[1]
  fields::image.plot(x = 1:50, y = 1:nStats, z = newMat, axes = FALSE, zlim = range, 
                     col = colorRampPalette(xCol)(30), xlab = "", ylab = "")
  
  abline(v = c(0,10,20,30,40,50) + 0.5)
  abline(v = c(0,10,20,30,40,50) - 1.5)
  abline(h = c(1:(nStats +1)) + 0.5)
  axis(2, at = 1:nrow(mat), labels = x.labels, las = 2, cex.axis = 1)
  axis(3, at = seq(5, 45, by = 10), labels = y.labels, cex.axis = 1.25, las = 1)
  axis(1, at = xpos, labels = rep(dimnames(mat)[[3]], 5), cex.axis = 0.9, las = 2)
  
  # text(x = xpos, y = 0, srt = 45, adj = 1, xpd = TRUE, labels = rep(dimnames(mat)[[3]], 5), cex = 0.65)
  box() 
  
  for(i in 1:dim(mat2)[1]){
    for(k in 1:dim(mat2)[2]){
       for(j in 1:dim(mat2)[3])
        if(!is.na(mat2[i,k,j]) & mat2[i,k,j] < 0.05){
          text(x = 10*(k-1) +j,  y = i, labels = "*" )
      }
    }
  }
}
```

```{r}
par(mar = c(5,10,5,3))

# Tidier labels
alt.x.labels = dimnames(R2MatFull)[[1]]

#NOTE: The following line only works if statisticsIndices was originally assigned statisticsIndices1 (with 19 elements)
alt.x.labels[c(7, 11, 14, 15, 16, 17, 18)] = c("Yule-PDA ratio", "mean I'", "MGL prin. eigenvalue", "MGL asymmetry", "MGL peakedness", "MGL eigengap", "nLTT stat")

# Modifying a few model abbreviations here (e.g. am and pw)
dimnames(R2MatFull)[[3]] = c("ca", "fh", "gen", "hs", "am", "pw", "ve", "xe")

image.real3d(R2MatFull, signMatFull, x.labels = alt.x.labels, y.labels = c("env. filtering", "dispersal", "niche conserv.", "mut./speciation", "competition"))

mtext("Models", 1, line = 3.5, cex = 2)
mtext("Metrics", 2, at = 10, cex = 2, line = 8)
mtext("Processes", 3, cex = 2, line = 3)

```

## Plot averaged across models

Define plotting function for 2d

IS THIS SECTION STILL NECESSARY? CURRENTLY THROWING ERRORS

```{r}

image.real <- function(mat, xCol = c("blue", "white", "white", "red"), 
                       range = c(-1,1), x.labels = rownames(mat), y.labels = colnames(mat)) { 
  mat <- t(mat)[,nrow(mat):1]
  fields::image.plot(mat, axes = FALSE, zlim = range, 
                     col = colorRampPalette(xCol)(30))
  axis(1, at = seq(0, 1, length = nrow(mat)), labels = x.labels)
  axis(2, at = seq(0, 1, length = ncol(mat)), labels = y.labels, las = 2)
  box() 
}


# 2-panel figure
par(mfrow = c(1, 2), mar = c(3,11,3,3))

# Values: % agreement, Color: average sign
image.real(signMat, x.labels = c('env', 'dis', 'nic', 'mut', 'com')) 

for(i in 1:length(statisticsIndices)){
  for(j in 1:length(predictorsIndices)){
    text((j-1)/(length(predictorsIndices)-1),
         1-(i-1)/(length(statisticsIndices)-1), 
         labels = 100*round(disAg2[i,j], digits = 2))
  }
}
mtext(paste0(rep("(", 5), numModelsPerPredictor, rep(")", 5)), 1, at = seq(0, 1, length = ncol(signMat)), line = 2)
title(main = "% agreement in effect direction")

# Values: average correlation coefficient, Color: average sign
image.real(signMat, x.labels = c('env', 'dis', 'nic', 'mut', 'com')) 
for(i in 1:length(statisticsIndices)){
  for(j in 1:length(predictorsIndices)){
    text((j-1)/(length(predictorsIndices)-1),
         1-(i-1)/(length(statisticsIndices)-1), 
         labels = round(R2Mat[i,j], digits = 2))
  }
}
mtext(paste0(rep("(", 5), numModelsPerPredictor, rep(")", 5)), 1, at = seq(0, 1, length = ncol(signMat)), line = 2)
title(main = "Average correlation coefficient")


image.real(R2Mat, range = c(-1,1), xCol = c("beige", "firebrick2")) 
title(main = "R2 between parameter and tree metric")

for(i in 1:length(statisticsIndices)){
  for(j in 1:length(predictorsIndices)){
    text((j-1)/(length(predictorsIndices)-1),
         1-(i-1)/(length(statisticsIndices)-1), 
         labels = round(R2Mat[i,j], digits = 2))
  }
}


image.real(disAg2) 
for(i in 1:length(statisticsIndices)){
  for(j in 1:length(predictorsIndices)){
    text((j-1)/(length(predictorsIndices)-1),
         1-(i-1)/(length(statisticsIndices)-1), 
         labels = round(disAg[i,j], digits = 2))
  }
}

```

# Model inversion

```{r}
library(ranger)
```

## Calculating RF and CV RF values

To be fair, we run this analysis using only the tree metrics that are not strongly correlated with richness

```{r, cache=cacheData}

# Tree statistics not strongly correlated with richness
statistics2 = colnames(simuData)[statisticsIndices2]

predictability2 = array(dim = c(length(statistics2), 
                                length(predictors), 
                                length(models)))

predictability2R2 = array(dim = c(length(predictors), 
                                  length(models)))
rownames(predictability2R2) = predictors
colnames(predictability2R2) = models
predictabilityR2CV = predictability2R2

for(j in 1:length(predictorsIndices)){
  for(k in 1:length(models)){
    tmp = simuData[simuData$model == models[k],c(predictorsIndices[j], statisticsIndices2)]
    tmp = tmp[complete.cases(tmp),]
    if(nrow(tmp) > 0){
      form = as.formula(paste(predictors[j], "~ ."))
      x = ranger(form, data = tmp, 
                 importance = 'permutation',
                scale.permutation.importance = T)
      predictability2[,j,k] = importance(x)     
      predictability2R2[j,k] = x$r.squared
      
      cvindices = sample(1:5, size = nrow(tmp), replace = T)
      r2 = rep(NA, 5)
      for(i in 1:5){
         x = ranger(form, data = tmp[cvindices != i, ])
         xp = predict(x, data = tmp[cvindices == i,] )
         r2[i] = cor(tmp[cvindices == i, 1], xp$predictions)
      }
      predictabilityR2CV[j,k] = mean(r2)
    } 
  }
}
```

## Predictability within a model

Internally calculated R2 values from the ranger package

```{r}
image.real(predictability2R2, range = c(0,1), xCol = c("beige", "firebrick2"))
title(main = "Random forest [R2] (non-richness metrics)")
```

The internally caluclated R2 values of a RF can sometimes be misleading. That's why we additionaly calculated R2 values from a 5-fold CV in the calculation loop above. Results are very similar

```{r}
image.real(predictabilityR2CV, range = c(0,1), xCol = c("beige", "firebrick2"))
title(main = "Random forest [R2] CV predictability (non-richness metrics)")
```

The same as a barplot

```{r}
par(mar = c(4, 6, 0, 0), mgp = c(4, 1, 0))

x = barplot(t(predictabilityR2CV), beside = T, las = 1, horiz = F, ylim = c(0,1.1), ylab = "Predicted-observed correlation", xaxt = "n", cex.lab = 1.8, cex.axis = 1.6)

x = as.vector(x)
x2 = as.vector(t(predictabilityR2CV))

axis(1, at = x[c(4, 12, 20, 28, 36)], labels = c("env.\nfiltering", "dispersal", "niche\nconserv.", "mutation/\nspeciation", "competition"), cex.axis = 1.6, line =1, tick = F)

for(i in 1:length(x)){
  if(!is.na(x2[i])){
    text(x[i],x2[i] + 0.02, modelAbbrevs$new[(i-1)%%8 + 1], srt = 90,cex= 1.4, adj = 0)    
  }
}

```

```{r}
predMean2 = apply(predictability2, c(1,2), mean, na.rm = T)
predSD2 = apply(predictability2, c(1,2), sd, na.rm = T)

rownames(predMean2) = statistics2
colnames(predMean2) = predictors

image.real(predMean2, range = NULL, xCol = c("beige", "firebrick2"))
title(main = "mean variable importance (non-richness metrics)")

for(i in 1:length(statisticsIndices2)){
  for(j in 1:length(predictorsIndices)){
    text((j-1)/(length(predictorsIndices)-1),
         1-(i-1)/(length(statisticsIndices2)-1), 
         labels = round(predMean2[i,j] / predSD2[i,j], digits = 2))
  }
}
```

## Between models

As a next step, we calculate predictability between models, i.e. we use the random forests tuned for one model to predict to the next model.

Cross-model predictability is measured by Kendall's rank correlation coefficient.

For the diagonal (within-model predictability), we use the CV predictability calculated in the previous step.

```{r}
predictors[2]

modelsTMP = models[-5]

CpredictabilityR2 = array(dim = c(length(modelsTMP), length(modelsTMP)))

for(k in 1:length(modelsTMP)){
  train = simuData[simuData$model == modelsTMP[k],c(predictorsIndices[2], statisticsIndices2)]
  train = train[complete.cases(train),]
  if(nrow(train) > 0){
    form = as.formula(paste(predictors[2], "~ ."))
    x = ranger(form, data = train)
  } 
  for(j in 1:length(modelsTMP)){
    test = simuData[simuData$model == modelsTMP[j],c(predictorsIndices[2], statisticsIndices2)]
    test = test[complete.cases(test),]
    if(nrow(test) > 0 & k != j){
      y = predict(x, test)
      CpredictabilityR2[k,j] = cor(y$predictions, test$dis1, method = "kendall")
    }
  }
}

# could add row and column sums 

rownames(CpredictabilityR2) = modelAbbrevs$new[modelAbbrevs$original %in% modelsTMP]
colnames(CpredictabilityR2) = modelAbbrevs$new[modelAbbrevs$original %in% modelsTMP]

tmp = predictabilityR2CV[2,-5]

for(i in 1:length(modelsTMP)) CpredictabilityR2[i,i] = tmp[i]

xCol = c("darkblue","blue", "lightblue", "gray94", "pink","red", "darkred")

par(mar = c(8,8,8,8))

image.real(CpredictabilityR2, range = c(-1,1), xCol = xCol)
title(main = paste("Cross-model predictability [r], mean =", round(mean(CpredictabilityR2), digits = 3)),xlab = "Model predicted to", ylab = "Model trained on\n", cex.lab = 1.5)

for(i in 1:length(modelsTMP)){
  for(j in 1:length(modelsTMP)){
    text((j-1)/(length(modelsTMP)-1),
         1-(i-1)/(length(modelsTMP)-1), 
         labels = round(CpredictabilityR2[i,j], digits = 2))
  }
}


#above
axis(3, labels = round(colMeans(abs(CpredictabilityR2)), digits = 2), 
     at = seq(0,1,length.out = 7), las = 2)

#right
axis(4, labels = round(rowMeans(abs(CpredictabilityR2)), digits = 2), 
     at = seq(0,1,length.out = 7), las = 2)

```

# Case studies

## within each model

Run random forests within each simulation model in order to predict model parameter values from simulation-derived tree shape metrics. Then plug in empirical tree shape metrics to infer a process parameter value assuming that the tree arose from the processes represented by a particular simulation model.

```{r}

# run a version of metricsForManyTrees() to get tree metric output for empirical sister clades (but that function assumes diff model ID naming convention, so had to manually tweak)
empiricalMetrics = read.table('empiricalSisterClades_treeOutput.txt', header = T, sep = '\t')
empiricalMetrics$clade = rep(c(1,2), 6)
empiricalMetrics$cladeName = c('Xenarthra', 'Afrotheria',
                               'Carnivora/\nPerissodactyla', 'Cetartiodactyla',
                               'Apodi', 'Trochili',
                               'Rodentia', 'Lagomorpha',
                               'Coraciiformes', 'Passeriformes',
                               'Pelecaniformes', 'Ciconiiformes')
empMetrics = empiricalMetrics[, 4:22]
empMetrics2 = empiricalMetrics[, c(7, 8, 12:22)]

predictedValues2 = data.frame(tree = character(),
                              model = character(), 
                              predictor = character(),
                              value = numeric())

predictability2 = array(dim = c(length(statistics2), length(predictors), length(models)))
predictability2R2 = array(dim = c(length(predictors), length(models)))

# For each combination of process (predictorsIndices) and simulation model, infer the strength of process for each empirical tree based on the RF

for(j in 1:length(predictorsIndices)){
  for(k in 1:length(models)){
    tmp = simuData[simuData$model == models[k],c(predictorsIndices[j], statisticsIndices2)]
    tmp = tmp[complete.cases(tmp),]
    if(nrow(tmp) > 0){
      form = as.formula(paste(predictors[j], "~ ."))
      x = ranger(form, data = tmp, importance = 'permutation',
                 scale.permutation.importance = T)
      predictions = predict(x, data = empMetrics2)
      
      tmpPredictions = data.frame(tree = word(empiricalMetrics$tree, 1, sep = fixed('.')),
                                  model = rep(models[k], nrow(empMetrics2)),
                                  predictor = rep(predictors[j], nrow(empMetrics2)),
                                  value = predictions$predictions)
      
      predictedValues2 = rbind(predictedValues2, tmpPredictions)
      
    } 
  }
}



# Need to standardize the predicted parameter values relative to the range of values examined within a given model.  
#   scaled = (x - min) / (max - min)          

#so first need to identify min and max values for each model-parameter combo:

modelParamRange = simuData %>%
  group_by(model) %>%
  summarize(env1Min = min(env1, na.rm = T),
            env1Max = max(env1, na.rm = T),
            com1Min = min(com1, na.rm = T),
            com1Max = max(com1, na.rm = T),
            dis1Min = min(dis1, na.rm = T),
            dis1Max = max(dis1, na.rm = T),
            nic1Min = min(nic1, na.rm = T),
            nic1Max = max(nic1, na.rm = T),
            mut1Min = min(mut1, na.rm = T),
            mut1Max = max(mut1, na.rm = T))

# Note that for some model/parameter combinations, the modeler varied parameter values uniformly on a log scale, while for others values varied uniformly on an arithmetic scale.

# I examined the skewness values of all sets of values, and determined that each of the following combinations had a skewness of 0.7 or greater and merited being log-transformed prior to scaling. 
# (Otherwise, the scaled values, e.g. in model 'hs' for dis1 or mut1 are close to 0)

# ca: mut1, com1
# hs: dis1, mut1
# xe: dis1, mut1

modelParamRange$com1Min[modelParamRange$model == 'ca'] = log10(modelParamRange$com1Min[modelParamRange$model == 'ca'])
modelParamRange$com1Max[modelParamRange$model == 'ca'] = log10(modelParamRange$com1Max[modelParamRange$model == 'ca'])
modelParamRange$dis1Min[modelParamRange$model %in% c('hs', 'xe')] = log10(modelParamRange$dis1Min[modelParamRange$model %in% c('hs', 'xe')])
modelParamRange$dis1Max[modelParamRange$model %in% c('hs', 'xe')] = log10(modelParamRange$dis1Max[modelParamRange$model %in% c('hs', 'xe')])
modelParamRange$mut1Min[modelParamRange$model %in% c('hs', 'xe', 'ca')] = log10(modelParamRange$mut1Min[modelParamRange$model %in% c('hs', 'xe', 'ca')])
modelParamRange$mut1Max[modelParamRange$model %in% c('hs', 'xe', 'ca')] = log10(modelParamRange$mut1Max[modelParamRange$model %in% c('hs', 'xe', 'ca')])

# Rearrange
modelParamMinMax = data.frame(model = rep(modelParamRange$model, 5), 
                              predictor = rep(c('env1', 'com1', 'dis1', 'nic1', 'mut1'), each = 8),
                              minVal = c(modelParamRange$env1Min, modelParamRange$com1Min, modelParamRange$dis1Min, 
                                         modelParamRange$nic1Min, modelParamRange$mut1Min),
                              maxVal = c(modelParamRange$env1Max, modelParamRange$com1Max, modelParamRange$dis1Max, 
                                         modelParamRange$nic1Max, modelParamRange$mut1Max))


# Also need to log-transform the same parameter values in predictedValues2
predictedValues2$value[predictedValues2$model == 'ca' & predictedValues2$predictor %in% c('mut1', 'com1')] = 
  log10(predictedValues2$value[predictedValues2$model == 'ca' & predictedValues2$predictor %in% c('mut1', 'com1')])

predictedValues2$value[predictedValues2$model %in% c('hs', 'xe') & predictedValues2$predictor %in% c('mut1', 'dis1')] = 
  log10(predictedValues2$value[predictedValues2$model %in% c('hs', 'xe') & predictedValues2$predictor %in% c('mut1', 'dis1')])


# Function for rescaling a value between the min and max such that the scaled value varies between 0 and 1.
# If either the min or max are NA, Inf, or -Inf, returns NA
scaleFunction = function(value, min, max) {
  
  if (!min %in% c(-Inf, NA, Inf) & !max %in% c(-Inf, NA, Inf)) {
    
    scaledValue = (value - min) / (max - min)
  
  } else {
  
    scaledValue = NA
  
  }
  return(scaledValue)
}




scaledPredictions2 = predictedValues2 %>%
  left_join(modelParamMinMax, by = c('model', 'predictor')) %>%
  mutate(pair = str_extract(tree, "[1-9]"),
         clade = rep(c(1, 2), 174))

scaledPredictions2$scaledValue = apply(scaledPredictions2[, c('value', 'minVal', 'maxVal')], 1, 
                             function(x) scaleFunction(value = x[1], min = x[2], max = x[3]))

#Okabe-Ito color blind friendly palette
modelColors = data.frame(model = unique(predictedValues2$model),
                         color = c(rgb(0,0,0),
                                   rgb(230/255, 159/255, 0),
                                   rgb(86/255, 180/255, 233/255),
                                   rgb(0, 158/255, 115/255),
                                   rgb(240/255, 228/255, 66/255),
                                   rgb(0, 114/255, 178/255),
                                   rgb(213/255, 94/255, 0),
                                   rgb(204/255, 121/255, 167/255)))

newModelLabels = data.frame(model = c('ca', 'fh', 'gen', 'hs', 'la', 'pontarp', 've', 'xe'), 
newLabel = c('ca', 'fh', 'gen', 'hs', 'am', 'pw', 've', 'xe'))

scaledPredictions2 = left_join(scaledPredictions2, modelColors, by = 'model') %>% 
  left_join(newModelLabels, by = 'model')



# Plot of model inversion for dispersal across 6 empirical sister clade pairs
# --bottom panel illustrates relative magnitude of mean.Iprime which was identified as the tree metric most diagnostic of dispersal

# Determines order of panels:
pairs = c(1, 2, 5, 4, 6, 8)

# Specify dispersal process
pro = 'dis1' 

par(mfrow = c(2, 3), mar = c(5, 3, 1, 0), oma = c(0, 3, 0, 0))

panelLetters = letters[1:6]
panelIndex = 0

for (p in pairs) {
  
  panelIndex = panelIndex + 1
    
  tmp = filter(scaledPredictions2, pair == p, predictor == pro)
    
  plot(c(0.5,2.5), range(tmp$scaledValue), type = 'n', xaxt = 'n', xlab = "", ylab = "", las = 1,
       ylim = c(0, 1), yaxt = 'n')
  
  mtext(empiricalMetrics$cladeName[empiricalMetrics$pair == p], 1, at = 1:2, line = 2.5, cex = 1.1, padj = 0.5)
  offset = .19
  mtext(c("mean I'", "VRD", "mean I'", "VRD"), 1, at = c(1.1 - offset, 1.05 + offset, 2.05 - offset, 2. + offset), line = .5, cex = .8)
  
  axis(2, at = seq(0.3, 1, by = 0.1), las = 1)
 
  abline(h = 0.25)
     
  for (m in unique(tmp$model)) {
      
    points(tmp$clade[tmp$model==m], tmp$scaledValue[tmp$model==m], type = 'b', col = tmp$color[tmp$model==m], lwd = 3)
    text(2.1, tmp$scaledValue[tmp$model==m & tmp$clade == 2], tmp$newLabel[tmp$model==m], adj = 0)
  }
  
  # panel label
  text(.55, 0.97, paste0("(", panelLetters[panelIndex], ")"), cex = 2)

  # Add barplots showing empirical values of mean.Iprime and VRD for the two sister clades
  par(new = T)
  
  tmpBarData = empiricalMetrics[empiricalMetrics$pair == p, c('mean.Iprime', 'VRD')] %>%
    rowwise() %>%
    mutate(quantile.Iprime = sum(mean.Iprime > simuData$mean.Iprime, na.rm = T)/nrow(simuData),
           quantile.VRD = sum(VRD > simuData$VRD, na.rm = T)/nrow(simuData)) %>%
    dplyr::select(quantile.Iprime, quantile.VRD) %>%
    t() %>%
    as.matrix()
  
  
  barplot(tmpBarData, beside = T, yaxt = 'n', ylim = c(0, 3.5), border = NA, xlim = c(0.25, 6.5), xaxt = 'n', space = c(0.05, 1))
  
  #barplot(as.matrix(tmpBarData[, 3:4]), 
  #        yaxt = 'n', ylim = c(0, 2.5), space = .95, border = NA, xlim = c(0.25, 4.5))  
  
  if (panelIndex %in% c(1, 4)) {
    mtext("Tree Metrics\n(quantiles)", 2, cex = .9, line = 1.8, at = 0.37)
    mtext("Inferred Dispersal", 2, cex = 1.1, line = 3, at = 2.2)
  }
  
}


```

# Supplemental Figures

```{r}

# Switching out 2 model abbreviations
simuData$model3 = simuData$model
simuData$model3 = gsub("pontarp", "pw", simuData$model3)
simuData$model3 = gsub("la", "am", simuData$model3)

par(mfrow = c(5, 4), mar = c(2, 2, 2, 1)); for (i in c(25,27:43)) { boxplot(simuData[,i] ~ simuData$model3, main = names(simuData)[i], col = rainbow(8))}
```

# Reproducibility information

```{r}
sessionInfo()
```