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.
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")
# It can be helpful to clear the working environment
rm(list = ls())
# Load the package
library(borealis)
tps
file formatcreate.tps(
input.filename = "wings.csv",
output.filename = "wings.tps",
id.factors = c('species','caste','digitizer'),
include.scale = TRUE,
invert.scale = TRUE)
tps
data into Rshapes <- read.tps("wings.tps")
Or use sample data available in the package
data("Bombus.forewings", package = "borealis")
shapes <- Bombus.forewings
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)
shapes <- align.reflect(shapes, top.pt = 2, left.pt = 1, links = fw.links )
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)
semiLMs <- matrix(
c(5,6,7, 14,15,16),
ncol = 3, byrow = TRUE
)
wing.gpa <- align.procrustes(shapes, curves = semiLMs)
wing.gpa <- align.procrustes(shapes, outlier.analysis = TRUE)
wing.gpa <- listed.gdf(wing.gpa)
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)
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)
write.provenance(wing.gpa,
output.filename = "~/Desktop/wing.shape.provenance.md",
title = "Preliminary wing shape data provenance")
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)
# 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")
# 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)
# 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)
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 <- gm.prcomp(size.model$residuals)
ggGMMplot(allometry.corrected.pca, group = wing.gpa$gdf$species,
group.title = 'species', convex.hulls = TRUE,
include.legend = TRUE)
# 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)
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)
# 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)