---
title: "Swallowtail Morphology"
author: "Hannah Owens"
date: "16 December, 2019"
output:
  pdf_document:
    toc: yes
  html_document:
    code_folding: hide
    df_print: paged
    number_sections: yes
    toc: yes
    toc_float: yes
  word_document:
    toc: yes
---

```{r setup, include = FALSE}
library(geomorph);
library(Evomorph);
library(colorRamps);
library(caper);
library(geiger);
library(phytools);
library(phylolm);
library(rr2);
library(MuMIn);
knitr::opts_chunk$set(echo = TRUE);
```

# Reading in the data

## Generating list of species of interest

First I read in the phylogeny, from Owens, *et al.* 2017.

```{r read in and plot the tree, results = "hide"}
phy <- read.tree("Appendix 1_DataForAnalysis/PapilioPhylogeny.nwk");
sppList <- phy$tip.label;
```

Next, gotta read in and clean those coordinates.

## Reading in Forewing Data

First, the forewing data get read in, and species represented by fewer than 10 specimens are removed.

```{r read in forewing data, results = "hide"}
coordsF <- readland.tps(file = "Appendix 1_DataForAnalysis/Forewing.tps", 
                        specID = "ID", 
                        readcurves = F);
#Remove scale bar, known crummy landmark, species that are not found in the phylogeny
coordsF <- coordsF[c(-1,-12:-13),,dimnames(coordsF)[[3]] %in% sppList] 
#Remove species with < 10 specimens
coordsF <- coordsF[,,!(dimnames(coordsF)[[3]] %in% 
                         names(table(dimnames(coordsF)[[3]])[table(dimnames(coordsF)[[3]])
                                                             <10]))] 
```

The forewing data is then aligned.

```{r align forewing data, results = "hide"}
#Landmark data Procrustes fit
alignedF <- gpagen(coordsF, print.progress = F);
summary(alignedF);
```

Next, outlier specimens are examined, and extreme outliers are removed (since these are likely erroneous).

```{r get rid of forewing outliers}
#Plot outliers
junk <- plotOutliers(alignedF$coords)#Ordered list of indices of specimens by "outlierness"

#Remove outliers (if necessary, i.e. red points
coordsF <- coordsF[,,-junk[1:2]]
alignedF <- gpagen(coordsF[,,], print.progress = F);
rm(coordsF, junk);
```
Only chucking one *Chilasa agestor* and one *Chilasa paradoxa* at this point.

Finally, here is the aligned forewing data.

```{r check out the resulting forewing alignment}
#Plot resulting data
plot(alignedF);
```

And here is what a plot of all the forewing data look like, color-coded by species.

```{r forewings in tangent space, fig.width = 7.5, fig.height = 7.5, results="hide"}
#Plotting tangent space
gp <- as.factor(dimnames(alignedF$coords)[[3]]); #Gets species IDs from landmark file
col.gp <- primary.colors(length(levels(gp)), 
                         steps = 5, 
                         no.white = TRUE)#Gets colors to assign to species
names(col.gp) <- levels(gp); #Assigns colors to names
col.gp <- col.gp[match(gp, names(col.gp))]; #Assigns colors to individual points
papTangentFAll <- plotTangentSpace(alignedF$coords, 
                                label=NULL, 
                                groups= col.gp, 
                                legend = F); #Does the actual plotting
rm(col.gp, gp)
```

## Reading in Hindwing Data

First, the hindwing data get read in, and species represented by fewer than 10 specimens are removed.

```{r read in hindwing data and remove species not represented in phylogeny, results = "hide"}
coordsH <- readland.tps(file = "Appendix 1_DataForAnalysis/Hindwing.tps", specID = "ID", 
                        readcurves = F);
#Remove scale bar, species that are not found in the phylogeny
coordsH <- coordsH[-13:-14,,dimnames(coordsH)[[3]] %in% unique(names(alignedF$Csize))] 
#Remove species represented by fewer than 10 specimens
coordsH <- coordsH[,,!(dimnames(coordsH)[[3]] %in%
                         names(table(dimnames(coordsH)[[3]])[table(dimnames(coordsH)[[3]])
                                                             <10]))] 
```

The hindwing data is then aligned.

```{r clean and align hindwing data, results = "hide"}
#Landmark data Procrustes fit
alignedH <- gpagen(coordsH[,,], print.progress = F);
summary(alignedH);
rm(coordsH);
```

Next, outlier specimens are examined, and extreme outliers are removed (since these are likely erroneous).

```{r get rid of hindwing outliers, results = "hide"}
#Plot outliers
junk <- plotOutliers(alignedH$coords)

#Remove outliers
alignedH <- gpagen(alignedH$coords[,,-junk[1]], print.progress = F);
```
One *Chilasa paradoxa* is a particularly distinct outlier and will be pitched.

Finally, here is the aligned hindwing data.

```{r check out the resulting hindwing alignment}
#Plot resulting data
plot(alignedH);
rm(junk)
```

And here is what a plot of all the hindwing data look like, color-coded by species.

```{r hindwings in tangent space, fig.width = 10, fig.height = 10, results="hide"}
#Plotting tangent space
gp <- as.factor(dimnames(alignedH$coords)[[3]]); #Gets species IDs from landmark file
col.gp <- primary.colors(length(levels(gp)), steps = 5, no.white = TRUE)
names(col.gp) <- levels(gp); #Assigns colors to names
col.gp <- col.gp[match(gp, names(col.gp))]; #Assigns colors to individual points
papTangentHAll <- plotTangentSpace(alignedH$coords, 
                                label=NULL, 
                                groups= col.gp, 
                                legend = F); #Does the actual plotting
rm(col.gp, gp);
```

## Checking allometries

Before we get any farther, it would be good to check to make sure shape variance is not being complicated by centroid size.

### Forewing

```{r forewing allometry, results = "asis"}
alignedFdFrame <- geomorph.data.frame(alignedF, species = dimnames(alignedF$coords)[[3]])
temp <- procD.lm(f1 = coords ~ log(Csize)*species, data = alignedFdFrame,
                 print.progress = F);
summary(temp)
gengp = c(24,24,21,24)
plot(temp, type = "regression", predictor = log(alignedFdFrame$Csize), 
     reg.type = "PredLine", 
     pch = gengp[as.numeric(as.factor(substr(dimnames(alignedFdFrame$coords)[[3]],1,3)))],
     bg = viridis::viridis(length(unique(alignedFdFrame$species)))[as.numeric(
       as.factor(alignedFdFrame$species))], xlab = "log(CS)",
     ylim = c(-.075,.15))
legend(9.5,.15, legend = c("Heraclides", "Non-Heraclides"), pch = c(21,24), bty = "o", 
       box.col = "gray")
rm(temp, alignedFdFrame)
```

The relationship between forewing size and shape is different in each species (and species has a much stronger effect on shape than size), so we we cannot easily separate the two.

### Hindwing
```{r hindwing allometry, results = "asis"}
alignedHdFrame <- geomorph.data.frame(alignedH, species = dimnames(alignedH$coords)[[3]]);
temp <- procD.lm(f1 = coords ~ log(Csize)*species, data = alignedHdFrame,
                 print.progress = F);
summary(temp)
gengp = c(24,24,21,24)
plot(temp, type = "regression", predictor = log(alignedHdFrame$Csize), 
     reg.type = "PredLine", 
     pch = gengp[as.numeric(as.factor(substr(dimnames(alignedHdFrame$coords)[[3]],1,3)))], 
     bg = viridis::viridis(length(unique(dimnames(alignedHdFrame$coords)[[3]])))[
       as.numeric(as.factor(dimnames(alignedHdFrame$coords)[[3]]))], 
     xlab = "log(CS)")
legend(0.0295,-.12, legend = c("Heraclides", "Non-Heraclides"), pch = c(21,24), bty = "o", 
       box.col = "gray")
rm(temp, alignedHdFrame)
```

As with forewings, the relationship between forewing size and shape is different in each species (and species has a much stronger effect on shape than size), so we we cannot easily separate the two.

## Calculating Mean Shape on a Per-Species Basis

Next, I calculated mean shapes for every species, as these are more appropriate for some of the comparative phylogenetic analyses.

### Forewing mean shape

First, I calculate the mean forewing shape for each species, and count how many specimens went into each mean shape. I also verify that at this stage all species are still represented by at least 10 specimens, since some specimens were removed as statistical outliers in earlier dataset cleaning steps.

```{r calculate mean forewing shape, results = "hide"}
sppList <- unique(dimnames(alignedF$coords)[[3]])
meanAlignF <- array(dim = c(10,2,length(sppList)));
measurementCount <- NULL;
for(i in 1:length(sppList)){
  meanAlignF[,,i] <- mshape(alignedF$coords[,,dimnames(alignedF$coords)[[3]]==sppList[[i]]]);
  measurementCount <- append(measurementCount,
                             length(alignedF$coords[,,dimnames(alignedF$coords)[[3]]==
                                                      sppList[[i]]])/20);
}
dimnames(meanAlignF)[[3]] <- sppList;
meanF <- gpagen(meanAlignF[,,measurementCount>=10], print.progress = F);
rm(measurementCount, meanAlignF, i)
```

Now I plot the outliers.

```{r get rid of mean forewing shape outliers, results = "hide"}
#Output will give you ordered list of indices of specimens
junk <- plotOutliers(meanF$coords)

#Remove outliers (if necessary, i.e. red points
meanAlignF <- meanF$coords[,,-junk[1]];
meanF <- gpagen(meanAlignF, print.progress = F);
rm(meanAlignF, junk)
```

Ok. Only took out *Pterourus euterpinus*.

And here is a plot of the results, color-coded by genus.

```{r plotting mean forewing shape, fig.width = 7.5, fig.height = 7.5, results = "hide"}
#Plotting tangent space
meanFphy <- drop.tip(phy, phy$tip[!phy$tip.label%in%dimnames(meanF$coords)[[3]]]);

# Assigning colors
col.gp <- c( "#5f3d98", "#fcb864", "#e56225", "#b2abd2");

#pdf("../Figures/ForewingPhyloMorphoSpace.pdf")
plotGMPhyloMorphoSpace(phy = meanFphy, A = meanF$coords, 
                       tip.labels = F, node.labels = F, 
                       ancStates = F, shadow = T,
                       plot.param = list(t.bg=c(rep(col.gp[1], times = 16), 
                                                rep(col.gp[2], times = 6), 
                                                rep(col.gp[3], times = 1),
                                                rep(col.gp[4], times = 25)),
                                         n.cex = 1, n.bg = c(rep("black", 3),
                                                        rep(unique(col.gp)[1], 15),
                                                        rep(unique(col.gp)[2], 5),
                                                        rep(unique(col.gp)[4], 24))));
legend(-0.067, -0.02, legend= c("Pterourus", "Chilasa", "Alexanoria", "Heraclides"), pch=19,
       col=col.gp, bty = "n")
#dev.off()
```

#### Mean shape allometry

Now it probably is worth doing another allometry test for the mean forewing dataset. 

```{r mean forewing allometry, results = "asis"}
meanFforAllo <- geomorph.data.frame(meanF, genus = substr(dimnames(meanF$coords)[[3]],1,3));
temp <- procD.lm(f1 = coords ~ log(Csize)*genus, data = meanFforAllo, print.progress = F);
summary(temp)
col.gp <- c("#e56225", "#fcb864", "#b2abd2", "#5f3d98");
plot(temp, type = "regression", predictor = log(meanFforAllo$Csize), reg.type = "PredLine", 
     pch = 21, bg = col.gp[as.numeric(as.factor(meanFforAllo$genus))], xlab = "log(CS)")
legend(0.0001, -0.028, legend= c("Alexanoria", "Chilasa", "Heraclides", "Pterourus"), pch=19,
       col=unique(col.gp), bty = "n")
rm(temp, meanFforAllo)
```
The relationship between forewing size and shape is pretty different in each genus (and genus has a much stronger effect on shape than size), so we we cannot easily separate the two.

### Hindwing Mean Shape

Next, I calculate the mean hindwing shape for each species, and count how many specimens went into each mean shape. I also verify that at this stage all species are still represented by at least 10 specimens, since some specimens were removed as statistical outliers in earlier dataset cleaning steps.

```{r calculate mean hindwing shape, results = "hide"}
meanAlignH <- array(dim = c(12,2,length(sppList)));
dimnames(meanAlignH)[[3]] <- sppList;
measurementCount <- NULL;
for(i in 1:length(sppList)){
  meanAlignH[,,i] <- mshape(alignedH$coords[,,dimnames(alignedH$coords)[[3]]== sppList[[i]]]);
  measurementCount <- append(measurementCount,
                             length(alignedH$coords[,,dimnames(alignedH$coords)[[3]]==
                                                      sppList[[i]]])/24);
}
meanH <- gpagen(meanAlignH[,,measurementCount>=10], print.progress = F);
rm(measurementCount, meanAlignH, i)
```


```{r get rid of mean hindwing shape outliers, results = "hide"}
junk <- plotOutliers(meanH$coords, inspect.outliers = T)
rm(junk);
```

There are no outliers.

And here is a plot of the results, color-coded by genus.

```{r plotting mean hindwing shape, fig.width = 7.5, fig.height = 7, results = "hide"}
meanHphy <- drop.tip(phy, phy$tip[!phy$tip.label%in%dimnames(meanH$coords)[[3]]]);

col.gp <- c( "#5f3d98", "#fcb864", "#e56225", "#b2abd2");

#pdf("../Figures/HindwingPhyloMorphoSpace.pdf")
plotGMPhyloMorphoSpace(phy = meanHphy, A = meanH$coords, tip.labels = F, 
                       node.labels = F, ancStates = F,
                       plot.param = list(t.bg=c(rep(col.gp[1], times = 16),
                                                rep(col.gp[2], times = 6),
                                                rep(col.gp[3], times = 1),
                                                rep(col.gp[4], times = 24)), 
                                         n.cex = 1, n.bg = c(rep("black", 3), 
                                                             rep(col.gp[1], 15), 
                                                             rep(col.gp[2], 5), 
                                                             rep(col.gp[4], 23))));
legend(-0.175, -0.1, legend= c("Pterourus", "Chilasa", "Alexanoria", "Heraclides"), 
       pch=19, col=unique(col.gp), bty = "n")
#dev.off()
rm(col.gp);
```

#### Mean shape allometry

Now it probably is worth doing another allometry test for the mean hindwing dataset. I guess.

```{r mean hindwing allometry, results = "asis"}
meanHforAllo <- geomorph.data.frame(meanH, 
                                    genus = substr(dimnames(meanH$coords)[[3]],1,3));
temp <- procD.lm(f1 = coords ~ log(Csize)*genus, data = meanHforAllo, print.progress = F);
summary(temp)
col.gp <- c("#e56225", "#fcb864", "#b2abd2", "#5f3d98");
plot(temp, type = "regression", predictor = log(meanHforAllo$Csize), reg.type = "PredLine",
     pch = 21, bg = col.gp[as.numeric(as.factor(meanHforAllo$genus))], xlab = "log(CS)")
legend(0, -0.11, legend= c("Alexanoria", "Chilasa", "Heraclides", "Pterourus"), pch=19,
       col=unique(col.gp), bty = "n")
rm(temp, meanHforAllo)
```

The relationship between forewing size and shape is different in each genus (and genus has a stronger effect on shape than size), so we cannot easily separate the two.

## Calculating other intraspecific wing shape metrics

We need four additional metrics each for fore- and hindwings at an intraspecific level. They are the mean procustes distance from the intraspecific mean shape, morphospace volume, mean centroid size, and centroid size variance. Morphospace volume is calculated by multiplying the ranges of each principal component (excluding any PC values more than 1.5 times greater or less than the interquartile range).  All four metrics are then centered and standardized as relative scores from 0 to 1.

### Forewing intraspecific shape metrics

First, we have to define a function that will rescale the data to values between 0 and 1, so that they are comparable when we are generating linear models.

```{r 0 to 1 min/max scaling function, echo=F}
rsc <- function (x) 
{
  nv <- vector(mode = "list", length = length(x));
  for (i in 1:length(x)){
    nv[[i]] <- (x[[i]] - min(x))/(max(x)-min(x))
  }
  return(unlist(nv));
}
```

Now we calculate the intraspecific shape metrics. Only the first seven principle components are used, which account for 95% of the variance in the dataset.

```{r calculating intraspecific forewing shape metrics, echo=FALSE, results='asis'}
#Species-level forewing shape metrics
cleanSpNames <- sort(unique(dimnames(alignedF$coords)[[3]]));
meanProcDist <- list()
meanCSize <- list()
varCSize <- list()
morphVol <- list()
sampleSize <- list()
for (i in 1:length(cleanSpNames)){
  temp0 <- alignedF$coords[,,dimnames(alignedF$coords)[[3]]==cleanSpNames[[i]]]
  if(is.null(dimnames(temp0))){
    meanProcDist[[i]] <- NA;
    meanCSize[[i]] <- NA;
    varCSize[[i]] <- NA;
    morphVol[[i]] <- NA;
    sampleSize[[i]] <- 1;
  }
  else{
    temp1 <- gpagen(temp0, print.progress = F)
    meanProcDist[[i]] <- mean(ShapeDist(temp1$coords,temp1$consensus))
    meanCSize[[i]] <- mean(temp1$Csize)
    varCSize[[i]] <- var(temp1$Csize)
    morphVol[[i]] <- prod(
      apply(papTangentFAll$pc.scores[rownames(papTangentFAll$pc.scores)==cleanSpNames[[i]],1:7], 
            FUN = function (x) quantile(x, probs=.95), MARGIN = 2) - 
        apply(papTangentFAll$pc.scores[rownames(papTangentFAll$pc.scores)==cleanSpNames[[i]],1:7], 
              FUN = function (x) quantile(x, probs=.05), MARGIN = 2)
      );
    sampleSize[[i]] <- length(dimnames(temp0)[[3]])
  }
}
forewingIntShMetrics <- as.data.frame(cbind(rsc(unlist(meanProcDist)), rsc(unlist(morphVol)),
                                            rsc(unlist(meanCSize)), rsc(unlist(varCSize)),
                                            unlist(sampleSize)));
row.names(forewingIntShMetrics) <- unlist(cleanSpNames);
colnames(forewingIntShMetrics) <- c("MeanProcrustesDist.", "Morph.Volume",
                                    "MeanCentroidSize", "CentroidSizeVar.", "SampleSize");
knitr::kable(forewingIntShMetrics, caption = "Intraspecific shape metrics for the forewing.");
rm(temp0, temp1, i, meanProcDist, morphVol, meanCSize, varCSize, sampleSize, papTangentFAll)
```

### Hindwing intraspecific shape metrics

As it says. Only the first eight principle components are used, which account for 95% of the variance in the dataset.

```{r calculating intraspecific hindwing shape metrics, echo=FALSE, results='asis'}
#Species-level hindwing shape metrics
cleanSpNames <- sort(unique(dimnames(alignedH$coords)[[3]]));
meanProcDist <- list()
meanCSize <- list()
varCSize <- list()
morphVol <- list()
sampleSize <- list()
for (i in 1:length(cleanSpNames)){
  temp0 <- alignedH$coords[,,dimnames(alignedH$coords)[[3]]==cleanSpNames[[i]]]
  temp1 <- gpagen(temp0, print.progress = F)
  meanProcDist[[i]] <- mean(ShapeDist(temp1$coords,temp1$consensus))
  meanCSize[[i]] <- mean(temp1$Csize)
  varCSize[[i]] <- var(temp1$Csize)
  morphVol[[i]] <- prod(
    apply(papTangentHAll$pc.scores[rownames(papTangentHAll$pc.scores)==cleanSpNames[[i]],1:7], 
          FUN = function (x) quantile(x, probs=.95), MARGIN = 2) - 
      apply(papTangentHAll$pc.scores[rownames(papTangentHAll$pc.scores)==cleanSpNames[[i]],1:7], 
            FUN = function (x) quantile(x, probs=.05), MARGIN = 2)
    );
  sampleSize[[i]] <- length(dimnames(temp0)[[3]])
}
hindwingIntShMetrics <- as.data.frame(cbind(rsc(unlist(meanProcDist)), rsc(unlist(morphVol)),
                                            rsc(unlist(meanCSize)), rsc(unlist(varCSize)),
                                            unlist(sampleSize)));
row.names(hindwingIntShMetrics) <- unlist(cleanSpNames);
colnames(hindwingIntShMetrics) <- c("MeanProcrustesDist.", "Morph.Volume",
                                    "MeanCentroidSize", "CentroidSizeVar.", "SampleSize");
knitr::kable(hindwingIntShMetrics, caption = "Intraspecific shape metrics for the hindwing.");
rm(temp0, temp1, i, meanProcDist, morphVol, meanCSize, varCSize, sampleSize, papTangentHAll)
```

# Analyses

Now that we have nice clean datasets, it's time to roll up our sleeves and crank through some analyses!

## Mean Shape

The first set of analyses are done on a per-species basis.

### Phylogenetic signal

Is there phylogenetic signal for forewing shape? More particularly, we are examining the hypothsesis that wing shape is evolutionarily constrained. This is predicated on the very simple idea that the more closely related two species are, the more similar they will look. To answer this question, we use Blomberg's K, which is the ratio of variance among species over the contrasts variance.

```{r phylo signal mean forewing shape, results = "asis"}
meanFphy <- drop.tip(phy, phy$tip[!phy$tip.label%in%dimnames(meanF$coords)[[3]]]);
meanFshapeSig <- physignal(A = meanF$coords, phy = meanFphy, iter = 10000, 
                           print.progress = F);
cat("Blomberg's K = ", round(meanFshapeSig$phy.signal, digits = 4), "; p = ", round(meanFshapeSig$pvalue, digits = 4), ".", sep = "");
```

For mean forewing shape, there is a fair amount of phylogenetic signal. 

Ok, but is there phylogenetic signal for hindwing shape?

```{r phylo signal mean hindwing shape, results = "asis"}
meanHphy <- drop.tip(phy, phy$tip[!phy$tip.label%in%dimnames(meanH$coords)[[3]]]);
meanHshapeSig <- physignal(A = meanH$coords, phy = meanHphy, iter = 10000, 
                           print.progress = F);
cat("Blomberg's K = ", round(meanHshapeSig$phy.signal, digits = 4), "; p = ", round(meanHshapeSig$pvalue, digits = 4), ".", sep = "");
```

For mean hindwing shape, there is moderate phylogenetic signal. Given these findings, it's suggestive that hindwing and forewing shapes are not constrained to evolve completely in tandem, but let's do a more direct test.

### Comparing evolutionary patterns in hindwing and forewing

The first step is merging the datasets in order to make the comparison.

```{r setting up mean fore- and hindwing evolutionary patterns, results = "hide"}
meanFforMerge <- meanF$coords[,,dimnames(meanF$coords)[[3]]%in%
                                dimnames(meanH$coords)[[3]]]
meanHforMerge <- meanH$coords[,,dimnames(meanH$coords)[[3]]%in%
                                dimnames(meanH$coords)[[3]]]
allTree <- drop.tip(meanFphy, meanFphy$tip[!meanFphy$tip.label%in%
                                             dimnames(meanHforMerge)[[3]]]);
combinedMeanMorphoData <- combine.subsets(meanFforMerge, meanHforMerge)
dimnames(combinedMeanMorphoData$coords)[[3]] <- dimnames(meanHforMerge)[[3]]
fhGrp <- c("f", "f", "f", "f", "f", "f", "f", "f", "f", "f", 
           "h", "h", "h", "h", "h", "h", "h", "h", "h", "h", "h", "h");
rm(meanFforMerge, meanHforMerge)
```

Then we can compare evolutionary rates of shapes across the tree.

```{r comparing mean evolutionary rate, results = "asis"}
#Evolutionary rate
compMF <- compare.multi.evol.rates(A = combinedMeanMorphoData$coords, gp = fhGrp, 
                                   Subset = F, phy = allTree, iter = 10000, 
                                   print.progress = F); 
cat("The evolutionary rate for the hindwing is: ", 
      round(compMF$sigma.d.gp[2], digits = 5), "; ", 
      "The evolutionary rate for the forewing is: ", 
      round(compMF$sigma.d.gp[1], digits = 5), "; ", 
      "The ratio between the two is: ", 
      round(compMF$sigma.d.ratio, digits = 4), "; ", 
      "P value: ", round(compMF$P.value, digits = 4), sep = "");
```



```{r testing phylogenetic integration, results = "asis"}
#Phylogenetic integration
phyloIntResults <- phylo.integration(combinedMeanMorphoData$coords[1:10,,], A2=combinedMeanMorphoData$coords[11:22,,], allTree, print.progress = F);
cat("rPLS: ", 
    round(phyloIntResults$r.pls, digits = 4), "; ", 
    "P value: ", round(phyloIntResults$P.value, 4), sep = "");
```

Based on this information, there is a fairly strong correlation between the variance in hindwing and forewing, corrected for the phylogeny, suggesting there is some balance required between hindwing and forewing shape. Based on the findings of these hypothesis tests, it seems that forewing evolution is fairly constrained compared to hindwing evolution.

***

## Csize of mean shape

Next, we will take a look at how forewing and hindwing centroid size evolve across the phylogeny, testing the same hypotheses we had for wing shape. This may suggest more about biomechanical balance in the system (or not!).

### Phylogenetic signal

As always, we check phylogenetic signal first. 

Here's the forewing.

```{r phylo signal mean forewing centroid size, results = "asis"}
meanFCsizeSig <- physignal(A = meanF$Csize, phy = meanFphy, iter = 10000, 
                           print.progress = F);
cat("Blomberg's K = ", round(meanFCsizeSig$phy.signal, digits = 4), "; p = ",
    round(meanFCsizeSig$pvalue, digits = 4), ".", sep = "");
```

And the hindwing.

```{r phylo signal mean hindwing centroid size, results = "asis"}
meanHCsizeSig <- physignal(A = meanH$Csize, phy = meanHphy, iter = 10000, 
                           print.progress = F);
cat("Blomberg's K = ", round(meanHCsizeSig$phy.signal, digits = 4), "; p = ",
    round(meanHCsizeSig$pvalue, digits = 4), ".", sep = "");
```

Based on this information, there is some significant phylogenetic signal in the forewing, but the hindwing is borderline.

### Plotting Centroid Sizes in Phylomorphospace
Next, we need to define a function to plot the centroid size results in phylomorphospace. This is an adaptation of `plotGMPhyloMorphoSpace()` from `geomorph`, but without doing a PCA on the data before doing the ancestral state reconstructions and plotting the results.

```{r function defininition for phlomorphospaceACEnoPC, include=FALSE}
phylomorphospaceACEnoPC <- function (phy, A, tip.labels = TRUE, node.labels = TRUE, 
                                     ancStates = TRUE, xaxis = 1, yaxis = 2, 
                                     plot.param = list(), shadow = FALSE) 
{
  # Safety checks
  if (any(is.na(A)) == T) {
    stop("Data matrix contains missing values. Estimate these first (see 'estimate.missing').")
  }
  if (length(dim(A)) == 3) {
    if (is.null(dimnames(A)[[3]])) {
      stop("Data matrix does not include taxa names as dimnames for 3rd dimension.")
    }
    x <- two.d.array(A)
  }
  if (length(dim(A)) == 2) {
    if (is.null(dimnames(A)[[1]])) {
      stop("Data matrix does not include taxa names as dimnames for rows.")
    }
    x <- A
  }
  if (!inherits(phy, "phylo")) 
    stop("tree must be of class 'phylo.'")
  if (!is.binary.tree(phy)) 
    stop("tree is not fully bifurcating (consider 'multi2di' in ape.")
  N <- length(phy$tip.label)
  Nnode <- phy$Nnode
  if (N != dim(x)[1]) {
    stop("Number of taxa in data matrix and tree are not not equal.")
  }
  if (length(match(rownames(x), phy$tip.label)) != N) 
    stop("Data matrix missing some taxa present on the tree.")
  if (length(match(phy$tip.label, rownames(x))) != N) 
    stop("Tree missing some taxa in the data matrix.")
  
  # Estimate ancestral states
  x <- x[phy$tip.label, ]
  anc.states <- NULL
  for (i in 1:ncol(x)) {
    x1 <- x[, i]
    tmp <- vector()
    for (j in 1:Nnode + N) {
      a <- multi2di(root(phy, node = j, resolve.root = F))
      tmp[j - N] <- ace(x1, a, method = "pic")$ace[1]
    }
    anc.states <- cbind(anc.states, tmp)
  }
  colnames(anc.states) <- NULL
  row.names(anc.states) <- 1:length(tmp)
  all.data <- rbind(x, anc.states)
  
  # Plotting
  p.p <- plot.param
  if (is.null(p.p$t.bg)) 
    p.p$t.bg = "black"
  if (is.null(p.p$t.pch)) 
    p.p$t.pch = 21
  if (is.null(p.p$t.cex)) 
    p.p$t.cex = 2
  if (is.null(p.p$n.bg)) 
    p.p$n.bg = "white"
  if (is.null(p.p$n.pch)) 
    p.p$n.pch = 21
  if (is.null(p.p$n.cex)) 
    p.p$n.cex = 1.25
  if (is.null(p.p$l.col)) 
    p.p$l.col = "black"
  if (is.null(p.p$lwd)) 
    p.p$lwd = 3
  if (is.null(p.p$txt.adj)) 
    p.p$txt.adj = c(-0.1, -0.1)
  if (is.null(p.p$txt.col)) 
    p.p$txt.col = "black"
  if (is.null(p.p$txt.cex)) 
    p.p$txt.cex = 1
  if (is.null(p.p$n.txt.adj)) 
    p.p$n.txt.adj = c(-0.1, -0.1)
  if (is.null(p.p$n.txt.col)) 
    p.p$n.txt.col = "black"
  if (is.null(p.p$n.txt.cex)) 
    p.p$n.txt.cex = 0.6

  if (tip.labels == TRUE) {
    plot(all.data[, xaxis], all.data[, yaxis], type = "n", 
         xlim = c(min(all.data[,xaxis]), max(all.data[,xaxis])), 
         ylim = c(min(all.data[,yaxis]), max(all.data[,yaxis])), 
         xlab = paste(xaxis), 
         ylab = paste(yaxis))
  }
  if (tip.labels == FALSE) {
    plot(all.data[, xaxis], all.data[, yaxis], type = "n",
         xlim = c(min(all.data[,xaxis]), max(all.data[,xaxis])), 
         ylim = c(min(all.data[,yaxis]), max(all.data[,yaxis])),
         xlab = paste(xaxis), ylab = paste(yaxis))
  }
  for (i in 1:nrow(phy$edge)) {
    lines(all.data[(phy$edge[i, ]), xaxis], 
          all.data[(phy$edge[i,]), yaxis], 
          type = "l", col = p.p$l.col, lwd = p.p$lwd)
  }
  points(all.data[1:N, xaxis], all.data[1:N, yaxis], pch = p.p$t.pch, 
         bg = p.p$t.bg, cex = p.p$t.cex)
  points(all.data[(N + 1):nrow(all.data), xaxis], 
         all.data[(N + 1):nrow(all.data), yaxis], 
         pch = p.p$n.pch, bg = p.p$n.bg, 
         cex = p.p$n.cex)
  if (tip.labels == TRUE) {
    text(all.data[1:N, xaxis], all.data[1:N, yaxis], rownames(all.data)[1:N], 
         col = p.p$txt.col, cex = p.p$txt.cex, adj = p.p$txt.adj)
  }
  if (node.labels == TRUE) {
    text(all.data[(N + 1):nrow(all.data), xaxis], 
         all.data[(N + 1):nrow(all.data), yaxis], 
         rownames(all.data)[(N + 1):nrow(all.data)], 
         col = p.p$n.txt.col, cex = p.p$n.txt.cex, 
         adj = p.p$n.txt.adj)
  }
  
  if (ancStates == TRUE) {
    return(anc.states)
  }
}
```

Now we plot.

```{r plotting centroid size with phylogeny}
centroidData <- cbind(meanF$Csize[names(meanF$Csize)%in%allTree$tip.label], 
                      meanH$Csize[names(meanH$Csize)%in%allTree$tip.label]);
colnames(centroidData) <- c("Forewing", "Hindwing");

col.gp <- c( "#5f3d98", "#fcb864", "#e56225", "#b2abd2");

#pdf("../Figures/CentroidSizePhyloMorphoSpace.pdf")
phylomorphospaceACEnoPC(phy = allTree, A = centroidData, tip.labels = F,
                        node.labels = F, ancStates = T, xaxis = 2, yaxis = 1,
                        plot.param = list(t.bg=c(rep(col.gp[1], 16),
                                                 rep(col.gp[2], 6),
                                                 rep(col.gp[3], 1),
                                                 rep(col.gp[4], 24)), 
                                          n.cex = 1, 
                                          n.bg = c(rep("black", 3), 
                                                   rep(col.gp[1], 15),
                                                   rep(col.gp[2], 5),
                                                   rep(col.gp[4], 23))))
legend(1.020, 1.00125, legend= c("Pterourus", "Chilasa", "Alexanoria", "Heraclides"), 
       pch=19, col=col.gp, bty = "n")
#dev.off()
```

### Evolutionary rate

In order to do this, I'll first fit a series of evolutionary models to the data. Then we'll see which is most appropriate, and also what the estimated evolutionary rate (sigma squared) under that model is.

As always, forewing first.

```{r forewing optimizing evolutionary model and estimating evolutionary rate, results = "asis"}
#Evolutionary rate
meanFCsize <- meanF$Csize;
forewingEvolModelCsize <- matrix(data = NA, nrow = 5, ncol = 3)
colnames(forewingEvolModelCsize) <- c("AICc", "Akaike Weight","Sigma Squared");
rownames(forewingEvolModelCsize) <- c("BM", "OU", "EB", "trend", "white");
alphaOU <- NULL
alphaCI <- NULL
count <- 1;
while (count <= nrow(forewingEvolModelCsize)){
  temp <- fitContinuous(meanFphy, meanFCsize, 
                        model = rownames(forewingEvolModelCsize)[[count]]);
  forewingEvolModelCsize[count, 1] <- temp$opt$aicc
  forewingEvolModelCsize[count, 3] <- temp$opt$sigsq
  if(count == 2){
    alphaOU <- temp$opt$alpha
    if(!(0 >= temp$bnd[1,1] && 0 <= temp$bnd[1,2])){
      alphaCI <- "does not include zero."
    }
    else{
      alphaCI <- "includes 0. DANGER!!"
    }
  }
  count <- count + 1;
}
forewingEvolModelCsize[,2] <- Weights(forewingEvolModelCsize[,1])
knitr::kable(forewingEvolModelCsize, caption = "Summary of evolutionary model testing 
             and evolutionary rate estimation for Papilio forewing centroid size.");

cat("\n\nThe evolutionary model that fits the forewing data best is: ",
    rownames(forewingEvolModelCsize)[forewingEvolModelCsize[,1]==min(forewingEvolModelCsize[,1])],
    ". \n\nThe sigma squared (an estimate of evolutionary rate) for forewing centroid size is: ",
    forewingEvolModelCsize[forewingEvolModelCsize[,1]==min(forewingEvolModelCsize[,1]),2], ".", 
    sep = "");

cat("\n\nThe alpha parameter under the OU model is: ", alphaOU, 
    ". The confidence interval of alpha ", alphaCI, sep = "");
```

So there is evolution going on and it's following an OU model. What about hindwing centroid size?

```{r hindwing optimizing evolutionary model and estimating evolutionary rate, results = "asis"}
#Evolutionary rate
meanHCsize <- meanH$Csize;
hindwingEvolModelCsize <- matrix(data = NA, nrow = 5, ncol = 3)
colnames(hindwingEvolModelCsize) <- c("AICc", "Akaike Weight","Sigma Squared");
rownames(hindwingEvolModelCsize) <- c("BM", "OU", "EB", "trend", "white");
alphaOU <- NULL
count <- 1;
while (count <= nrow(hindwingEvolModelCsize)){
  temp <- fitContinuous(meanHphy, meanHCsize, 
                        model = rownames(hindwingEvolModelCsize)[[count]]);
  hindwingEvolModelCsize[count, 1] <- temp$opt$aicc
  hindwingEvolModelCsize[count, 3] <- temp$opt$sigsq
  if(count == 2){
    alphaOU <- temp$opt$alpha
    if(!(0 >= temp$bnd[1,1] && 0 <= temp$bnd[1,2])){
      alphaCI <- "does not include zero."
    }
    else{
      alphaCI <- "includes 0. DANGER!!"
    }
  }
  count <- count + 1;
}
hindwingEvolModelCsize[,2] <- Weights(hindwingEvolModelCsize[,1])
knitr::kable(hindwingEvolModelCsize, 
             caption = "Summary of evolutionary model testing and evolutionary rate 
             estimation for Papilio hindwing centroid size.");

cat("\n\nThe evolutionary model that fits the hindwing data best is: ",
    rownames(hindwingEvolModelCsize)[hindwingEvolModelCsize[,1]==min(hindwingEvolModelCsize[,1])],
    ". \n\nThe sigma squared (an estimate of evolutionary rate) for hindwing centroid size is: ",
    hindwingEvolModelCsize[hindwingEvolModelCsize[,1]==min(hindwingEvolModelCsize[,1]),2], ".", 
    sep = "");

cat("\n\nThe alpha parameter under the OU model is: ", alphaOU, 
    ". The confidence interval of alpha ", alphaCI, sep = "");
```

Well look at that! While white noise fits best, OU is very close. This is probably not surprising, since there was less phylogenetic signal in the hindwing than the forewing. It is also looks like the hindwing is evolving much faster than the forewing. This makes sense, since some groups have tails and some don't.

Finally, let's check the phylogenetic integration of the fore- and hindwing (note that this comparison is done assuming Brownian motion). This is a statistical test comparing the significance of the difference we note above.

```{r phylogenetic integration of forewing and hindwing centroid size}
meanFforMerge <- meanFCsize[names(meanFCsize) %in% names(meanHCsize)]
meanHforMerge <- meanHCsize[names(meanHCsize) %in% names(meanFforMerge)]
allCentroidSize <- data.frame(names(meanFforMerge), meanFforMerge,meanHforMerge)
colnames(allCentroidSize) <- c("names","forewingCsize", "hindwingCsize")
coCZVatMeanCsize <- comparative.data(allTree, data = allCentroidSize, names.col = "names");

CMEvolRatesCS <- compare.multi.evol.rates(coCZVatMeanCsize$data, coCZVatMeanCsize$phy, 
                                          Subset = F, gp = c("f", "h"), print.progress = F);

cat("\n\nThe evolutionary rate for the hindwing is: ", 
      round(CMEvolRatesCS$sigma.d.gp[2], digits = 6), "; ", 
      "\n\nThe evolutionary rate for the forewing is: ", 
      round(CMEvolRatesCS$sigma.d.gp[1], digits = 8), "; ", 
      "\n\nThe ratio between the two is: ", 
      round(CMEvolRatesCS$sigma.d.ratio, digits = 6), "; ", 
      "\n\nP value: ", round(CMEvolRatesCS$P.value, digits = 6), sep = "");
```

Yep. Hindwing is significantly faster.

### Correlation between fore- and hindwing centroid size

There's one last thing that may be worth trying--is there a correlation between fore- and hindwing centroid size?

First, I define a pair of functions to calculate the AICc scores for three alternative evolutionary models (Brownian motion, OU (fixed root), and early burst) and return the model with the best fit based on minimum AICc score.

```{r functions for phylolm model fit-finding}
modelTestPhyloLM <- function(formula, data, phy){
  phylolmModels <- c("BM", "OUfixedRoot", "EB");
  modScores <- vector(mode = "numeric", length(phylolmModels))
  for (i in 1:length(phylolmModels)){
    tmp <- phylolm(formula = formula, data = data, phy = phy, model = phylolmModels[i]);
    modScores[i] <- MuMIn::AICc(tmp)
  }
  names(modScores) <- phylolmModels
  return(modScores)
}

bfphylolm <- function(scores){
  return(names(sort(scores)[1]))
}
```

Now how about that correlation between fore- and hindwing centroid size?

```{r forewing and hindwing centroid size correlation, results="asis"}
#Correlation
phylolmFHCsizeEModels <- modelTestPhyloLM(formula = hindwingCsize ~ forewingCsize, 
                                          data = coCZVatMeanCsize$data, 
                                          phy = coCZVatMeanCsize$phy)
phylolmFHCsize <- phylolm(formula = hindwingCsize ~ forewingCsize, 
                          data = coCZVatMeanCsize$data, 
                          phy = coCZVatMeanCsize$phy, 
                          model = bfphylolm(phylolmFHCsizeEModels));
cat("\n\nThe correlation coefficient of hindwing size as explained by forewing size is: ",
    phylolmFHCsize$coefficients[2], ", *p* = ", summary(phylolmFHCsize)$coefficients[2,4], 
    ". The $R^2_pred$ value for the correlation between hindwing and forewing centroid size is ",
    R2.pred(phylolmFHCsize, phy = coCZVatMeanCsize$phy), ". ", 
    "\n\nThis is based on the best-fit evolutionary model, which was ",
    bfphylolm(phylolmFHCsizeEModels), ".", sep = "");
```

***

## Taxonomically-separate analyses

It is possible that for the two major clades of "New World" Papilio butterflies, the lowland tropical Heraclides and that highland tropical/temperate Pterourus + Agehana + Chilasa + Alexanoria, different selective forces may be at work. This is especially likely for the second, which is the group with many species that are mimics.

### Mean Shape

The first set of analyses are done on a per-species basis.

#### Phylogenetic signal

Is there phylogenetic signal for forewing shape? More particularly, we are examining the hypothsesis that wing shape is evolutionarily constrained. This is predicated on the very simple idea that the more closely related two species are, the more similar they will look. To answer this question, we use Blomberg's K, which is the ratio of variance among species over the contrasts variance.

```{r taxonomically split phylo signal mean forewing shape, results = "asis"}
meanFphyHer <- extract.clade(phy = meanFphy, node = 72);
meanFtrimHer <- gpagen(meanF$coords[,,dimnames(meanF$coords)[[3]]%in%
                                      meanFphyHer$tip.label],
                       print.progress = F);
meanFshapeSigHer <- physignal(A = meanFtrimHer$coords, phy = meanFphyHer, iter = 10000,
                              print.progress = F);

meanFphyPter <- extract.clade(phy = meanFphy, node = 50);
meanFtrimPter <- gpagen(meanF$coords[,,dimnames(meanF$coords)[[3]]%in%
                                       meanFphyPter$tip.label],
                        print.progress = F);
meanFshapeSigPter <- physignal(A = meanFtrimPter$coords, phy = meanFphyPter, iter = 10000,
                               print.progress = F);

#Plot comparisons of all signals
fwsignals <- c(round(meanFshapeSig$phy.signal, digits = 4), 
               round(meanFshapeSigHer$phy.signal, digits = 4), 
               round(meanFshapeSigPter$phy.signal, digits = 4))
fwpvals <- c(round(meanFshapeSig$pvalue, digits = 4), round(meanFshapeSigHer$pvalue, digits = 4),
             round(meanFshapeSigPter$pvalue, digits = 4))
compFWSigRes <- cbind(fwsignals, fwpvals)
rownames(compFWSigRes) <- c("All_Taxa", "Heraclides", "Non-Heraclides");
colnames(compFWSigRes) <- c("Blomberg's K", "P-Value");
knitr::kable(compFWSigRes, 
             caption = "Summary of phylogenetic signal of forewing shape, 
             broken down by clade.");
```

Looking at each of the taxa individually, *Heraclides* has a lower signal than the other two. Cool cool cool. What about for the hindwing?

```{r taxonomically split phylo signal mean hindwing shape, results = "asis"}
meanHphyHer <- extract.clade(phy = meanHphy, node = 71);
meanHtrimHer <- gpagen(meanH$coords[,,dimnames(meanH$coords)[[3]]%in%
                                      meanHphyHer$tip.label],
                       print.progress = F);
meanHshapeSigHer <- physignal(A = meanHtrimHer$coords, phy = meanHphyHer, iter = 10000,
                              print.progress = F);

meanHphyPter <- extract.clade(phy = meanHphy, node = 49);
meanHtrimPter <- gpagen(meanH$coords[,,dimnames(meanH$coords)[[3]]%in%
                                       meanHphyPter$tip.label],
                        print.progress = F);
meanHshapeSigPter <- physignal(A = meanHtrimPter$coords, phy = meanHphyPter, iter = 10000,
                               print.progress = F);

#Plot comparisons of all signals
hwsignals <- c(round(meanHshapeSig$phy.signal, digits = 4), 
               round(meanHshapeSigHer$phy.signal, digits = 4), 
               round(meanHshapeSigPter$phy.signal, digits = 4))
hwpvals <- c(round(meanHshapeSig$pvalue, digits = 4), 
             round(meanHshapeSigHer$pvalue, digits = 4),
             round(meanHshapeSigPter$pvalue, digits = 4))
compHWSigRes <- cbind(hwsignals, hwpvals)
rownames(compHWSigRes) <- c("All_Taxa", "Heraclides", "Non-Heraclides");
colnames(compHWSigRes) <- c("Blomberg's K", "P-Value");
knitr::kable(compHWSigRes, 
             caption = "Summary of phylogenetic signal of hindwing shape, 
             broken down by clade.");
```

For mean hindwing shape, there is moderate phylogenetic signal, but when the taxa are broken up, it's clear that there is much stronger phylogenetic signal for shape in non-Heraclides species than Heraclides species. This is weird, since the non-Heraclides were expected to have *less* signal because that's where all the mimics are. But, when you think of it in the context of Ackerley 2009 (PNAS), you can interpret this as showing less convergent evolution. This, I think, helps picture what's going on...

#### Comparing evolutionary patterns in hindwing and forewing

The first step is merging the datasets in order to make the comparison.

```{r setting up mean fore- and hindwing taxon split evolutionary patterns, results = "hide"}
#Heraclides
meanFforMergeHer <- gpagen(meanFtrimHer$coords[,,dimnames(meanFtrimHer$coords)[[3]]%in%
                                                 dimnames(meanHtrimHer$coords)[[3]]],
                           print.progress = F)
meanHforMergeHer <- gpagen(meanHtrimHer$coords[,,dimnames(meanHtrimHer$coords)[[3]]%in%
                                                 dimnames(meanFtrimHer$coords)[[3]]],
                           print.progress = F)
allTreeHer <- drop.tip(phy, phy$tip[!phy$tip.label%in%dimnames(meanHforMergeHer$coords)[[3]]]);
combinedMeanMorphoDataHer <- combine.subsets(meanFforMergeHer, meanHforMergeHer)
dimnames(combinedMeanMorphoDataHer$coords)[[3]] <- dimnames(meanHforMergeHer$coords)[[3]]

#Pterourus
meanFforMergePter <- gpagen(meanFtrimPter$coords[,,dimnames(meanFtrimPter$coords)[[3]]%in%
                                                   dimnames(meanHtrimPter$coords)[[3]]],
                            print.progress = F)
meanHforMergePter <- gpagen(meanHtrimPter$coords[,,dimnames(meanHtrimPter$coords)[[3]]%in%
                                                   dimnames(meanFtrimPter$coords)[[3]]],
                            print.progress = F)
allTreePter <- drop.tip(phy, phy$tip[!phy$tip.label%in%dimnames(meanHforMergePter$coords)[[3]]]);
combinedMeanMorphoDataPter <- combine.subsets(meanFforMergePter, meanHforMergePter)
dimnames(combinedMeanMorphoDataPter$coords)[[3]] <- dimnames(meanHforMergePter$coords)[[3]]
```

Then we can compare evolutionary rates of shapes across the tree.

```{r taxon split comparing mean evolutionary rate, results = "asis"}
#Evolutionary rate
compMFHer <- compare.multi.evol.rates(A = combinedMeanMorphoDataHer$coords, 
                                      gp = fhGrp, Subset = F, phy = allTreeHer, iter = 10000,
                                      print.progress = F);
compMFPter <- compare.multi.evol.rates(A = combinedMeanMorphoDataPter$coords, 
                                       gp = fhGrp, Subset = F, phy = allTreePter, iter = 10000,
                                       print.progress = F);

#Putting everything together for reporting
forewingRates <- c(round(compMF$sigma.d.gp[1], digits = 7), 
                   round(compMFHer$sigma.d.gp[1], digits = 7), 
                   round(compMFPter$sigma.d.gp[1], digits = 7));
hindwingRates <- c(round(compMF$sigma.d.gp[2], digits = 7), 
                   round(compMFHer$sigma.d.gp[2], digits = 7), 
                   round(compMFPter$sigma.d.gp[2], digits = 7));
ratioFtoH <- c(round(compMF$sigma.d.ratio, digits = 7),
               round(compMFHer$sigma.d.ratio, digits = 7),
               round(compMFPter$sigma.d.ratio, digits = 7));
pValues <- c(compMF$P.value,
             compMFHer$P.value,
             compMFPter$P.value);

compEvShapeTaxSplit <- cbind(forewingRates, hindwingRates, ratioFtoH, pValues);
colnames(compEvShapeTaxSplit) <- c("Forewing_Rate", "Hindwing_Rate", "Ratio_H_to_F", "p_Value");
rownames(compEvShapeTaxSplit) <- c("All_Taxa", "Heraclides", "Non-Heraclides");

knitr::kable(compEvShapeTaxSplit, 
             caption = "Summary of evolutionary rate comparisons between forewing 
             and hindwing, broken down by clade.");
```

Based on this information, it appears that in both groups, hindwings are evolving more quickly than hindwings. This is more pronounced in the Heraclides taxa than the non-Heraclides taxa. 

```{r testing phylogenetic integration with taxonomic split, results = "asis"}
#Phylogenetic integration
phyloIntResultsHer <- phylo.integration(combinedMeanMorphoDataHer$coords[1:10,,],
                                        A2=combinedMeanMorphoDataHer$coords[11:22,,], 
                                        allTreeHer, print.progress = F);
phyloIntResultsPter <- phylo.integration(combinedMeanMorphoDataPter$coords[1:10,,],
                                         A2=combinedMeanMorphoDataPter$coords[11:22,,],
                                         allTreePter, print.progress = F);

#Making the table
integrationResults <- c(round(phyloIntResults$r.pls, digits = 4),
                        round(phyloIntResultsHer$r.pls, digits = 4),
                        round(phyloIntResultsPter$r.pls, digits = 4));
integrationPvals <- c(round(phyloIntResults$P.value, digits = 4),
                      round(phyloIntResultsHer$P.value, digits = 4),
                      round(phyloIntResultsPter$P.value, digits = 4));
phyloIntTable <- cbind(integrationResults,integrationPvals);
rownames(phyloIntTable) <- c("All_Taxa", "Heraclides", "Non-Heraclides");
colnames(phyloIntTable) <- c("rPLS", "P value")

knitr::kable(phyloIntTable, caption = "Summary of phylogenetic integration test 
             between forewing and hindwing, broken down by clade.");
```

Things are very consistent. There is strong correlation between fore- and hindwing in all cases.

***

### Mean Csize

Next, we're back to looking at how centroid size evolves across the phylogeny, testing the same hypotheses we had for wing shape.

#### Phylogenetic signal

As always, we check phylogenetic signal first. 

Here's the forewing.

```{r phylo signal mean forewing centroid size with taxonomic split, results = "asis"}
#Heraclides
meanFCsizeSigHer <- physignal(A = meanFtrimHer$Csize, phy = meanFphyHer, iter = 10000,
                              print.progress = F);

#Pterourus
meanFCsizeSigPter <- physignal(A = meanFtrimPter$Csize, phy = meanFphyPter, iter = 10000,
                               print.progress = F);

#Summarize the results and print them
sigMeanFCsize <- c(round(meanFCsizeSig$phy.signal, digits = 4),
                       round(meanFCsizeSigHer$phy.signal, digits = 4),
                       round(meanFCsizeSigPter$phy.signal, digits = 4));
pValuesMeanFCsize <- c(round(meanFCsizeSig$pvalue, digits = 4),
                       round(meanFCsizeSigHer$pvalue, digits = 4),
                       round(meanFCsizeSigPter$pvalue, digits = 4));
summaryMeanFCsize <- cbind(sigMeanFCsize, pValuesMeanFCsize);
colnames(summaryMeanFCsize) <- c("Blomberg's K", "P value");
rownames(summaryMeanFCsize) <- c("All_Taxa", "Heraclides", "Non-Heraclides");

knitr::kable(summaryMeanFCsize, 
             caption = "Summary of phylogenetic signal of forewing centroid size, 
             broken down by clade.");
```

Interestingly, there is more and more statistically significant phylogenetic signal when the dataset is considered as a whole.

And the hindwing.

```{r phylo signal mean hindwing centroid size with taxonomic split, results = "asis"}
#Heraclides
meanHCsizeSigHer <- physignal(A = meanHtrimHer$Csize, phy = meanHphyHer, iter = 10000,
                              print.progress = F);

#Pterourus
meanHCsizeSigPter <- physignal(A = meanHtrimPter$Csize, phy = meanHphyPter, iter = 10000,
                               print.progress = F);

#Summarize the results and print them
sigmeanHCsize <- c(round(meanHCsizeSig$phy.signal, digits = 4),
                       round(meanHCsizeSigHer$phy.signal, digits = 4),
                       round(meanHCsizeSigPter$phy.signal, digits = 4));
pValuesmeanHCsize <- c(round(meanHCsizeSig$pvalue, digits = 4),
                       round(meanHCsizeSigHer$pvalue, digits = 4),
                       round(meanHCsizeSigPter$pvalue, digits = 4));
summarymeanHCsize <- cbind(sigmeanHCsize, pValuesmeanHCsize);
colnames(summarymeanHCsize) <- c("Blomberg's K", "P value");
rownames(summarymeanHCsize) <- c("All_Taxa", "Heraclides", "Non-Heraclides");

knitr::kable(summarymeanHCsize, 
             caption = "Summary of phylogenetic signal of hindwing centroid size, 
             broken down by clade.");
```

As with the forewing, the least centroid size significance is in *Heraclides*. However, in this case there's stronger, but less significant, signal in the non-*Heraclides* clade compared to all taxa.

#### Evolutionary rate

I wonder, are the best-fit models of evolutionary rate the same for each of the taxa, and how about the rates?

As always, forewing first.

```{r forewing optimizing evolutionary model and estimating evolutionary rate, taxonomically split, results = "asis"}
#Heraclides Evolutionary rate
forewingEvolModelCsizeHer <- matrix(data = NA, nrow = 5, ncol = 3)
colnames(forewingEvolModelCsizeHer) <- c("AICc", "Akaike Weight","Sigma Squared");
rownames(forewingEvolModelCsizeHer) <- c("BM", "OU", "EB", "trend", "white");
alphaOU <- NULL
count <- 1;
while (count <= nrow(forewingEvolModelCsizeHer)){
  temp <- fitContinuous(meanFphyHer, meanFtrimHer$Csize, 
                        model = rownames(forewingEvolModelCsizeHer)[[count]]);
  forewingEvolModelCsizeHer[count, 1] <- temp$opt$aicc
  forewingEvolModelCsizeHer[count, 3] <- temp$opt$sigsq
  if(count == 2){
    alphaOU <- temp$opt$alpha
    if(!(0 >= temp$bnd[1,1] && 0 <= temp$bnd[1,2])){
      alphaCI <- "does not include zero."
    }
    else{
      alphaCI <- "includes 0. DANGER!!"
    }
  }
  count <- count + 1;
}
forewingEvolModelCsizeHer[,2] <- Weights(forewingEvolModelCsizeHer[,1])
knitr::kable(forewingEvolModelCsizeHer, 
             caption = "Summary of evolutionary model testing and evolutionary rate 
             estimation for Heraclides forewing centroid size.");

cat("\n\nThe evolutionary model that fits the Heraclides forewing data best is: ", rownames(forewingEvolModelCsize)[forewingEvolModelCsize[,1]==min(forewingEvolModelCsize[,1])], 
    ". The sigma squared (an estimate of evolutionary rate) for forewing centroid size is: ",
    forewingEvolModelCsizeHer[forewingEvolModelCsizeHer[,1]==min(forewingEvolModelCsizeHer[,1]),2],
    ".\n", sep = "");
cat("\n\nThe alpha parameter under the OU model is: ", alphaOU, 
    ". The confidence interval of alpha ", alphaCI, sep = "");

#Pterourus Evolutionary rate
forewingEvolModelCsizePter <- matrix(data = NA, nrow = 5, ncol = 3)
colnames(forewingEvolModelCsizePter) <- c("AICc", "Akaike Weight","Sigma Squared");
rownames(forewingEvolModelCsizePter) <- c("BM", "OU", "EB", "trend", "white");
alphaOU <- NULL
count <- 1;
while (count <= nrow(forewingEvolModelCsizePter)){
  temp <- fitContinuous(meanFphyPter, meanFtrimPter$Csize, 
                        model = rownames(forewingEvolModelCsizePter)[[count]]);
  forewingEvolModelCsizePter[count, 1] <- temp$opt$aicc
  forewingEvolModelCsizePter[count, 3] <- temp$opt$sigsq
  if(count == 2){
    alphaOU <- temp$opt$alpha
    if(!(0 >= temp$bnd[1,1] && 0 <= temp$bnd[1,2])){
      alphaCI <- "does not include zero."
    }
    else{
      alphaCI <- "includes 0. DANGER!!"
    }
  }
  count <- count + 1;
}
forewingEvolModelCsizePter[,2] <- Weights(forewingEvolModelCsizePter[,1])
knitr::kable(forewingEvolModelCsizePter, 
             caption = "Summary of evolutionary model testing and evolutionary rate 
             estimation for Pterourus forewing centroid size.");

cat("\n\nThe evolutionary model that fits the Pterourus forewing data best is: ",
    rownames(forewingEvolModelCsize)[forewingEvolModelCsize[,1]==min(forewingEvolModelCsize[,1])],
    ". The sigma squared (an estimate of evolutionary rate) for forewing centroid size is: ",
    forewingEvolModelCsizePter[forewingEvolModelCsizePter[,1]==
                                 min(forewingEvolModelCsizePter[,1]),2], ".\n", sep = "");
cat("\n\nThe alpha parameter under the OU model is: ", alphaOU, 
    ". The confidence interval of alpha ", alphaCI, sep = "");
```

The best-fit evolutionary model for centroid size of the forewing for both Heraclides and Pterourus are consistent with the full clade, which is to say, it's all random.

```{r hindwing optimizing evolutionary model and estimating evolutionary rate, taxonomically split, results = "asis"}
#Heraclides Evolutionary rate
hindwingEvolModelCsizeHer <- matrix(data = NA, nrow = 5, ncol = 3)
colnames(hindwingEvolModelCsizeHer) <- c("AICc", "Akaike Weight","Sigma Squared");
rownames(hindwingEvolModelCsizeHer) <- c("BM", "OU", "EB", "trend", "white");
alphaOU <- NULL
count <- 1;
while (count <= nrow(hindwingEvolModelCsizeHer)){
  temp <- fitContinuous(meanHphyHer, meanHtrimHer$Csize, 
                        model = rownames(hindwingEvolModelCsizeHer)[[count]]);
  hindwingEvolModelCsizeHer[count, 1] <- temp$opt$aicc
  hindwingEvolModelCsizeHer[count, 3] <- temp$opt$sigsq
  if(count == 2){
    alphaOU <- temp$opt$alpha
    if(!(0 >= temp$bnd[1,1] && 0 <= temp$bnd[1,2])){
      alphaCI <- "does not include zero."
    }
    else{
      alphaCI <- "includes 0. DANGER!!"
    }
  }
  count <- count + 1;
}
hindwingEvolModelCsizeHer[,2] <- Weights(hindwingEvolModelCsizeHer[,1])
knitr::kable(hindwingEvolModelCsizeHer, 
             caption = "Summary of evolutionary model testing and evolutionary rate 
             estimation for Heraclides hindwing centroid size.");

cat("\n\nThe evolutionary model that fits the Heraclides hindwing data best is: ",
    rownames(hindwingEvolModelCsize)[hindwingEvolModelCsize[,1]==min(hindwingEvolModelCsize[,1])], 
    ". The sigma squared (an estimate of evolutionary rate) for hindwing centroid size is: ",
    hindwingEvolModelCsizeHer[hindwingEvolModelCsizeHer[,1]==min(hindwingEvolModelCsizeHer[,1]),2],
    ".\n", sep = "");

cat("\n\nThe alpha parameter under the OU model is: ", alphaOU, 
    ". The confidence interval of alpha ", alphaCI, sep = "");

#Pterourus Evolutionary rate
hindwingEvolModelCsizePter <- matrix(data = NA, nrow = 5, ncol = 3)
colnames(hindwingEvolModelCsizePter) <- c("AICc", "Akaike Weight","Sigma Squared");
rownames(hindwingEvolModelCsizePter) <- c("BM", "OU", "EB", "trend", "white");
alphaOU <- NULL
count <- 1;
while (count <= nrow(hindwingEvolModelCsizePter)){
  temp <- fitContinuous(meanHphyPter, meanHtrimPter$Csize, 
                        model = rownames(hindwingEvolModelCsizePter)[[count]]);
  hindwingEvolModelCsizePter[count, 1] <- temp$opt$aicc
  hindwingEvolModelCsizePter[count, 3] <- temp$opt$sigsq
  if(count == 2){
    alphaOU <- temp$opt$alpha
    if(!(0 >= temp$bnd[1,1] && 0 <= temp$bnd[1,2])){
      alphaCI <- "does not include zero."
    }
    else{
      alphaCI <- "includes 0. DANGER!!"
    }
  }
  count <- count + 1;
}
hindwingEvolModelCsizePter[,2] <- Weights(hindwingEvolModelCsizePter[,1])
knitr::kable(hindwingEvolModelCsizePter, 
             caption = "Summary of evolutionary model testing and evolutionary rate estimation 
             for Pterourus hindwing centroid size.");

cat("\n\nThe evolutionary model that fits the Pterourus hindwing data best is: ", rownames(hindwingEvolModelCsizePter)
    [hindwingEvolModelCsizePter[,1]==min(hindwingEvolModelCsizePter[,1])], 
    ". The sigma squared (an estimate of evolutionary rate) for hindwing centroid size is: ",
    hindwingEvolModelCsizePter[hindwingEvolModelCsizePter[,1]==
                                 min(hindwingEvolModelCsizePter[,1]),2], ".\n", sep = "");

cat("\n\nThe alpha parameter under the OU model is: ", alphaOU, 
    ". The confidence interval of alpha ", alphaCI, sep = "");
```

Yep, also nothing but randomness.

```{r integration of forewing and hindwing centroid size with clades split}
#Heraclides
meanFforMergeHer <- meanFtrimHer$Csize[names(meanFtrimHer$Csize) %in% names(meanHtrimHer$Csize)]
meanHforMergeHer <- meanHtrimHer$Csize[names(meanHtrimHer$Csize) %in% names(meanFforMergeHer)]
allCentroidSizeHer <- data.frame(names(meanFforMergeHer), meanFforMergeHer, meanHforMergeHer);
colnames(allCentroidSizeHer) <- c("names","forewingCsize", "hindwingCsize");
coCZVatMeanCsizeHer <- comparative.data(phy, data = allCentroidSizeHer, names.col = "names");
tempHer <- compare.multi.evol.rates(coCZVatMeanCsizeHer$data, coCZVatMeanCsizeHer$phy, 
                                    Subset = F, gp = c("f", "h"), print.progress = F);

#Pterourus
meanFforMergePter <- meanFtrimPter$Csize[names(meanFtrimPter$Csize) 
                                         %in% names(meanHtrimPter$Csize)]
meanHforMergePter <- meanHtrimPter$Csize[names(meanHtrimPter$Csize) 
                                         %in% names(meanFforMergePter)]
allCentroidSizePter <- data.frame(names(meanHforMergePter), meanFforMergePter, meanHforMergePter);
colnames(allCentroidSizePter) <- c("names","forewingCsize", "hindwingCsize");
coCZVatMeanCsizePter <- comparative.data(phy, data = allCentroidSizePter, names.col = "names");
tempPter <- compare.multi.evol.rates(coCZVatMeanCsizePter$data, coCZVatMeanCsizePter$phy, 
                                     Subset = F, gp = c("f", "h"), print.progress = F);
```

```{r Making the summary table for phylo integration}
forewingSigma <- c(round(CMEvolRatesCS$sigma.d.gp[1],9), 
                   round(tempHer$sigma.d.gp[1],9), 
                   round(tempPter$sigma.d.gp[1],9));
hindwingSigma <- c(CMEvolRatesCS$sigma.d.gp[2], tempHer$sigma.d.gp[2], tempPter$sigma.d.gp[2]);
ratio <- c(CMEvolRatesCS$sigma.d.ratio, tempHer$sigma.d.ratio, tempPter$sigma.d.ratio);
pValsPhyInt <- c(CMEvolRatesCS$P.value, tempHer$P.value, tempPter$P.value);
summaryEvolIntegration <- cbind(round(forewingSigma,9), 
                                hindwingSigma, ratio, pValsPhyInt);
colnames(summaryEvolIntegration) <- c("Forewing Rate", "Hindwing Rate", "Ratio", "P-value")
rownames(summaryEvolIntegration) <- c("All_Taxa", "Heraclides", "Non-Heraclides");

knitr::kable(summaryEvolIntegration, 
             caption = "Summary of phylogenetic integration tests for centroid size.");
```

Yep. Hindwing is significantly faster.

#### Correlation between fore- and hindwing centroid size

There's one last thing to test--what happens when you split out the two taxa and and look at the size correlation.

```{r forewing and hindwing centroid size correlation with taxon splits, results="asis"}
#Heraclides
phylolmFHCsizeHerEModels <- modelTestPhyloLM(formula = hindwingCsize ~ forewingCsize, 
                                             data = coCZVatMeanCsizeHer$data, 
                                             phy = coCZVatMeanCsizeHer$phy)
phylolmFHCsizeHer <- phylolm(formula = hindwingCsize ~ forewingCsize, 
                             data = coCZVatMeanCsizeHer$data, 
                             phy = coCZVatMeanCsizeHer$phy, 
                             model = bfphylolm(phylolmFHCsizeHerEModels));
cat("\n\nThe $R^2_pred$ value for the correlation between Heraclides hindwing and forewing 
    centroid size is ", 
    R2.pred(phylolmFHCsizeHer, phy = coCZVatMeanCsizeHer$phy), ". ", 
    "\n\nThis is based on the best-fit evolutionary model, which was ",
    bfphylolm(phylolmFHCsizeHerEModels), ".", sep = "");

#Pterourus
phylolmFHCsizePterEModels <- modelTestPhyloLM(formula = hindwingCsize ~ forewingCsize, 
                                              data = coCZVatMeanCsizePter$data, 
                                              phy = coCZVatMeanCsizePter$phy)
phylolmFHCsizePter <- phylolm(formula = hindwingCsize ~ forewingCsize, 
                              data = coCZVatMeanCsizePter$data, 
                              phy = coCZVatMeanCsizePter$phy, 
                              model = bfphylolm(phylolmFHCsizePterEModels));
cat("\n\nThe $R^2_pred$ value for the correlation between Pterourus hindwing and forewing 
    centroid size is ", 
    R2.pred(phylolmFHCsizePter, phy = coCZVatMeanCsizePter$phy), ". ", 
    "This is based on the best-fit evolutionary model, which was ",
    bfphylolm(phylolmFHCsizePterEModels), ".", sep = "");

#Printing out the results
cSizeCorrelation <- c(R2.pred(phylolmFHCsize, phy = coCZVatMeanCsize$phy), 
                      R2.pred(phylolmFHCsizeHer, phy = coCZVatMeanCsizeHer$phy),
                      R2.pred(phylolmFHCsizePter, phy = coCZVatMeanCsizePter$phy));
betas <- c(phylolmFHCsize$coefficients[2], 
           phylolmFHCsizeHer$coefficients[2], 
           phylolmFHCsizePter$coefficients[2]);
pvals <- c(summary(phylolmFHCsize)$coefficients[2,4],
           summary(phylolmFHCsizeHer)$coefficients[2,4],
           summary(phylolmFHCsizePter)$coefficients[2,4])

cSizeCorrelationTable <- cbind(cSizeCorrelation, betas, pvals);
colnames(cSizeCorrelationTable) <- c("Rsq Pred", "Forewing Beta", "P value for Beta");
rownames(cSizeCorrelationTable) <- c("All_Taxa", "Heraclides", "Non-Heraclides");

knitr::kable(cSizeCorrelationTable, 
             caption = "Summary of R^2 pred for forewing and hindwing centroid size with 
             taxa split out.");
```

Well there you go. Pretty consistent. Fairly low R squareds. It is weird that there was a stronger beta for all taxa than there was for either sub-clade, but only the overall clade had a significant coeffiecient. Therefore, there may just not have been enough samples in the subclades to detect the overall pattern.