knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.