This document is a quick starting point for running an analysis of geometric morphometrics using the borealis package. For a more detailed explanation of each step, please see the accompanying tutorial.

Installing the borealis package

# If you have a previous version, it's good to detach the old one first
detach("package:borealis", unload = TRUE)

devtools::install_github("aphanotus/borealis")

Initializing

# It can be helpful to clear the working environment
rm(list = ls())

# Load the package
library(borealis)

Convert raw landmark coordinates into the tps file format

create.tps(
  input.filename = "wings.csv",
  output.filename = "wings.tps",
  id.factors = c('species','caste','digitizer'),
  include.scale = TRUE,
  invert.scale = TRUE)

Import tps data into R

shapes <- read.tps("wings.tps")

Or use sample data available in the package

data("Bombus.forewings", package = "borealis")
shapes <- Bombus.forewings

View the shape data

fw.links <- matrix(c(1,2, 1,5, 5,4, 4,3, 3,2, 5,6, 6,7, 7,8, 8,9, 9,4, 3,11, 11,12, 11,10, 9,10, 10,14, 14,15, 15,16, 16,18, 18,20, 16,17, 17,8, 12,13, 13,19, 14,13, 18,19, 2,12),
                   ncol = 2, byrow = TRUE)

landmark.plot(shapes, links = fw.links)
landmark.plot(shapes, specimen.number = 1:4, links = fw.links)

Reflect specimens

shapes <- align.reflect(shapes, top.pt = 2, left.pt = 1, links = fw.links )

Angular alignment

data("mantis", package = "borealis")

# Define mantis.lines
{
  x <- 1:16
  mantis.lines <- matrix(c(x[-length(x)],x[-1]), ncol = 2)
  mantis.lines[10,] <- c(10,1)
  mantis.lines[15,] <- c(15,6)
  mantis.lines <- rbind(mantis.lines,
                        matrix(c(5,11, 6,11, 13,16, 14,16), ncol = 2, byrow = TRUE))
}

landmark.plot(mantis, specimen.number = 1:3, panels = c(1,3), links = mantis.lines)

mant.2 <- align.angle(mantis,
                      art.pt = 11,
                      angle.pts.1 = c(1:10),
                      angle.pts.2 = c(12:15),
                      rot.pts = c(12:16) )

landmark.plot(mant.2, specimen.number = 1:3, panels = c(1,3), links = mantis.lines)

Procrustes alignment

GPA with semilandmarks

semiLMs <- matrix(
  c(5,6,7,  14,15,16),
  ncol = 3, byrow = TRUE
)

wing.gpa <- align.procrustes(shapes, curves = semiLMs)

GPA with outlier detection

wing.gpa <- align.procrustes(shapes, outlier.analysis = TRUE)

The geomorph data frame

wing.gpa <- listed.gdf(wing.gpa)

Subsetting shape data

x <- which(wing.gpa$gdf$caste=="W")
workers <- subsetGMM(wing.gpa, specimens = x)

distal.wing <- subsetGMM(wing.gpa, landmarks = 1:10)
landmark.plot(distal.wing)

Custom additions to data provenance

wing.gpa$gdf$digitizer[1] <- "FJ"

wing.gpa <- add.provenance(
  wing.gpa,
  name="error.correction",
  title = "Corrected error in the metadata",
  text = "Upon checking lab notes, an error needed to be corrected. The digitizer of specimen 1 appears to have been FJ, not JL." )

cat(wing.gpa$provenance$error.correction)

Reporting data provenance

write.provenance(wing.gpa, 
                 output.filename = "~/Desktop/wing.shape.provenance.md",
                 title = "Preliminary wing shape data provenance")

Ordination

PCA

wing.pca <- gm.prcomp(wing.gpa$gdf$coords)
summary(wing.pca)

ggGMMplot(wing.pca, group = wing.gpa$gdf$species, 
          group.title = 'species', 
          convex.hulls = TRUE, include.legend = TRUE)

ggGMMplot(wing.pca, group = wing.gpa$gdf$species, 
          group.title = 'species', convex.hulls = TRUE,
          backtransform.examples = TRUE,
          ref.shape = mshape(wing.gpa$gdf$coords),
          shape.method = "TPS",
          bt.shape.mag = 3,
          bt.links = fw.links)

Phylogenetic PCA

# filter to just workers
x <- which(wing.gpa$gdf$caste=="W")
workers <- subsetGMM(wing.gpa, specimens = x)

# Find mean shape for each species
coords.by.species <- coords.subset(workers$gdf$coords, group = workers$gdf$species)
mshape.by.species <- lapply(coords.by.species, mshape)
species.names <- names(mshape.by.species)
mshape.by.species <- array(
  data = unlist(mshape.by.species),
  dim = c(dim(mshape.by.species[[1]])[1], 2, length(species.names)),
  dimnames = list(NULL,NULL,species.names)
)

# Prep the phylogeny
data("Bombus.tree", package = "borealis")
btree <- Bombus.tree

btree$tip.label <- btree$code.name

# Check that species.names (our abbreviations in the shape data) are all present in the tree
species.names %in% btree$tip.label

library(phytools)
btree <- keep.tip(btree, species.names)

# Simple overlay of the phylogeny onto the PC-space
pca.w.phylo <- gm.prcomp(mshape.by.species, phy = btree)
plot(pca.w.phylo, phylo = TRUE, main = "PCA with phylogeny")

# Phylogenetically-aligned PCA
paca <- gm.prcomp(mshape.by.species, phy = btree, align.to.phy = TRUE)
plot(paca, phylo = TRUE, main = "Phylogenetically-aligned PCA")

Modeling

# Set the number of iterations
# Higher values increase accuracy as well as compute time
i <- 1e3-1 

# A simple allometric model of shape
size.model <- procD.lm(coords ~ log(Csize), data = wing.gpa$gdf, iter = i) 
anova(size.model)

# A model with size and species
species.model <- procD.lm(coords ~ log(Csize) + species, data=wing.gpa$gdf, iter=i) 
anova(species.model) 

# Model comparison
anova(size.model, species.model)

Post hoc pairwise comparisons

# Filter to the 4 most common species, to ensure decent sample sizes
common.sp <- subsetGMM(wing.gpa, specimens = (wing.gpa$gdf$species %in% c("bimac","imp","tern","vag")))

# Build null and full models
common.size.model <- procD.lm(coords ~ log(Csize), data=common.sp$gdf, iter=i) 
anova(common.size.model)
common.sp.model <- procD.lm(coords ~ log(Csize) + species, data=common.sp$gdf, iter=i) 
anova(common.sp.model)

# Pairwise comparisons
common.sp.pw <- pairwise(fit = common.sp.model, 
                         fit.null = common.size.model,
                         groups = common.sp$gdf$species)
summary(common.sp.pw)

Testing for common allometry

species.unique.model <- procD.lm(coords ~ log(Csize) * species, data=wing.gpa$gdf, iter=i) 

# geomorph's allometry visualization
plotAllometry(fit = species.unique.model, size = log10(wing.gpa$gdf$Csize), 
              col = as.factor(wing.gpa$gdf$species), xlab = "log10 centroid size")

# the equivalent using borealis
gg.scaling.plot(
  x = log10(wing.gpa$gdf$Csize), y = wing.pca$x[,1],
  group=wing.gpa$gdf$species, group.title = "species",
  xlab = "log10 wing centroid size", ylab = "wing shape PC1",
  include.legend = TRUE,
  groups.trendlines = TRUE,
  fixed.aspect = FALSE
)

# Compare the models
anova(species.model, species.unique.model)

Allometry-corrected PCA

allometry.corrected.pca <- gm.prcomp(size.model$residuals)
ggGMMplot(allometry.corrected.pca, group = wing.gpa$gdf$species, 
          group.title = 'species', convex.hulls = TRUE,
          include.legend = TRUE)

Phylogenetic ANOVA

# Create a custom geomorph.data.frame with the relevant data
species.gdf <- geomorph.data.frame(
  coords = mshape.by.species,
  Csize = c(by(workers$gdf$Csize, workers$gdf$species, mean)),
  tree = btree
)

size.pgls <- procD.pgls(coords ~ log(Csize), phy = tree, data = species.gdf, iter = i)

anova(size.pgls)

Disparity comparisons

x <- which(workers$gdf$species %in% c("bimac","imp","tern"))
well.sampled.workers <- subsetGMM(workers, specimens = x)

morphol.disparity(coords ~ Csize, groups = ~ species, data = well.sampled.workers$gdf)

Modularity testing

# Define a modularity hypothesis
modularity.hypothesis1 <- rep("proximal",wing.gpa$landmark.number)
modularity.hypothesis1[c(1,4,5,6:9,17)] <- "distal"

# Test the hypothesis
mt1 <- modularity.test(wing.gpa$gdf$coords, modularity.hypothesis1, CI=TRUE, iter=99 )
summary(mt1)