--- title: "Analyses repeated with breeding season-specific foraging data." output: html_document --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE) ``` ## Introduction This document will run through all analyses described in our manuscript, only using the breeding-specific foraging information. Details of what is being done, rationale, and a walk-through of the results are not available in full here but can be found in the markdown documents "1_-Foraging_Ecology_and_Shape.Rmd/html" and "3--Variability_of_Allometry_Across_Families.Rmd/html". We have included a few brief annotations to indicate which model results are reported in our supplementary information and to help indicate when a new analysis beginning. ```{r message=FALSE} load("Charadriiform_GeoMorph_Data_breeding.RData") library(Morpho) library(dplyr) library(geomorph) library(motmot) library(ggplot2) library(ape) ``` ## 1:The influence of guild and size on charadriiform skull shape We will begin with the analysis on the whole skull. ```{r wholegdf} table(gdf.b$guild_narrow) singletons <- names(which(table(gdf.b$guild_narrow)==1)) singletons <- names(gdf.b$guild_narrow[match(singletons, gdf.b$guild_narrow)]) gdf.b_narrow<-gdf.b gdf.b_narrow$phy <- drop.tip(gdf.b_narrow$phy, singletons) gdf.b_narrow$Coords<-gdf.b_narrow$Coords[,,-match(singletons, dimnames(gdf.b_narrow$Coords)[[3]])] gdf.b_narrow$mass <-gdf.b_narrow$mass[-match(singletons, names(gdf.b_narrow$mass))] gdf.b_narrow$logCS <-gdf.b_narrow$logCS[-match(singletons, names(gdf.b_narrow$logCS))] gdf.b_narrow$Family <-gdf.b_narrow$Family[-match(singletons, names(gdf.b_narrow$Family))] gdf.b_narrow$guild_broad <-gdf.b_narrow$guild_broad[-match(singletons, names(gdf.b_narrow$guild_broad))] gdf.b_narrow$guild_narrow <-gdf.b_narrow$guild_narrow[-match(singletons, names(gdf.b_narrow$guild_narrow))] ``` ```{r wholeskullANOVA} narrowguild_model<- procD.lm(Coords~mass*guild_narrow, SS.type= "II", data= gdf.b_narrow, warpgrids = FALSE, iter=999, print.progress = FALSE) summary(narrowguild_model) residuals <-narrowguild_model$residuals residuals.pca<- prcomp(residuals) PCscores.residuals <-residuals.pca$x[,1:37] lambdaestimate <-transformPhylo.ML(y=PCscores.residuals, phy=gdf.b_narrow$phy, model= "lambda", lowerBound = 0, upperBound = 1) lambdaestimate$Lambda ``` The MLE of lamba is 0, so we will report the non-phylogenetic approach. If we want to be really conservative, we could take upper bound estimate for lambda, just to see what would happen: ```{r wholeskullupper} transformedtree<- geiger:::rescale(gdf.b_narrow$phy, "lambda", 0.0018) cformvcv <- vcv(transformedtree) narrowguild_model.phylo <- procD.lm(Coords~mass*guild_narrow, SS.type= "II",Cov=cformvcv, data= gdf.b_narrow, warpgrids = FALSE, iter=999, print.progress = FALSE) summary(narrowguild_model.phylo) summary(narrowguild_model) ``` Now, we will repeat these analyses on the beak and braincase seperately. ```{r seperatelms} braincase.lms <- c(1:4,17,20,21,32:35, 39:70) beak.lms <-c(1:255)[-braincase.lms] beak.coords <-gdf.b$Coords[c(beak.lms),,] dim(beak.coords) <- c(212, 3, 262) dimnames(beak.coords)[[3]] <- dimnames(gdf.b$Coords)[[3]] ``` Beginning with the beak. ```{r beakgdf} gdf.b.beak <-geomorph.data.frame(Coords=beak.coords, phy= gdf.b$phy, mass=gdf.b$mass, Family = gdf.b$Family, guild_broad= gdf.b$guild_broad, guild_narrow = gdf.b$guild_narrow) table(gdf.b.beak$guild_narrow) singletons <- names(which(table(gdf.b.beak$guild_narrow)==1)) singletons <- names(gdf.b.beak$guild_narrow[match(singletons, gdf.b.beak$guild_narrow)]) gdf.b_narrow.beak<-gdf.b.beak gdf.b_narrow.beak$phy <- drop.tip(gdf.b_narrow.beak$phy, singletons) gdf.b_narrow.beak$Coords<-gdf.b_narrow.beak$Coords[,,-match(singletons, dimnames(gdf.b_narrow.beak$Coords)[[3]])] gdf.b_narrow.beak$mass <-gdf.b_narrow.beak$mass[-match(singletons, names(gdf.b_narrow.beak$mass))] gdf.b_narrow.beak$Family <-gdf.b_narrow.beak$Family[-match(singletons, names(gdf.b_narrow.beak$Family))] gdf.b_narrow.beak$guild_broad <-gdf.b_narrow.beak$guild_broad[-match(singletons, names(gdf.b_narrow.beak$guild_broad))] gdf.b_narrow.beak$guild_narrow <-gdf.b_narrow.beak$guild_narrow[-match(singletons, names(gdf.b_narrow.beak$guild_narrow))] ``` ```{r beakANOVA} narrowguild_model.beak <- procD.lm(Coords~mass*guild_narrow, SS.type= "II", data= gdf.b_narrow.beak, warpgrids = FALSE, iter=999, print.progress = FALSE) summary(narrowguild_model.beak) ``` ```{r beaklambda} residuals.beak <-narrowguild_model.beak$residuals residuals.pca.beak<- prcomp(residuals.beak) PCscores.residuals.beak <-residuals.pca.beak$x[,1:37] lambdaestimate.beak <-transformPhylo.ML(y=PCscores.residuals.beak, phy=gdf.b_narrow.beak$phy, model= "lambda", lowerBound = 0, upperBound = 1) lambdaestimate.beak$Lambda ``` Again, this is about 0 so we report the non-phylogenetic version, but we can also try the upper estimate just to be safe. ```{r beakupper} transformedtree.beak<- geiger:::rescale(gdf.b_narrow.beak$phy, "lambda", lambdaestimate.beak$Lambda[,3]) cformvcv <- vcv(transformedtree.beak) narrowguild_model.beak_phylo<- procD.lm(Coords~mass*guild_narrow, SS.type= "II",Cov=cformvcv, data= gdf.b_narrow.beak, warpgrids = FALSE, iter=999, print.progress = FALSE) summary(narrowguild_model.beak_phylo) summary(narrowguild_model.beak) ``` Now, the braincase. ```{r brcgdf} braincase.coords<- gdf.b$Coords[c(braincase.lms),,] dim(braincase.coords) <- c(43, 3, 262) dimnames(braincase.coords)[[3]] <- dimnames(gdf.b$Coords)[[3]] #assign the names to the matrix of coordinates gdf.b.brc <-geomorph.data.frame(Coords=braincase.coords, phy= gdf.b$phy, mass=gdf.b$mass, Family = gdf.b$Family, guild_broad= gdf.b$guild_broad, guild_narrow = gdf.b$guild_narrow) singletons <- names(which(table(gdf.b.brc$guild_narrow)==1)) singletons <- names(gdf.b.brc$guild_narrow[match(singletons, gdf.b.brc$guild_narrow)]) gdf.b_narrow.brc<-gdf.b.brc gdf.b_narrow.brc$phy <- drop.tip(gdf.b_narrow.brc$phy, singletons) gdf.b_narrow.brc$Coords<-gdf.b_narrow.brc$Coords[,,-match(singletons, dimnames(gdf.b_narrow.brc$Coords)[[3]])] gdf.b_narrow.brc$mass <-gdf.b_narrow.brc$mass[-match(singletons, names(gdf.b_narrow.brc$mass))] gdf.b_narrow.brc$Family <-gdf.b_narrow.brc$Family[-match(singletons, names(gdf.b_narrow.brc$Family))] gdf.b_narrow.brc$guild_broad <-gdf.b_narrow.brc$guild_broad[-match(singletons, names(gdf.b_narrow.brc$guild_broad))] gdf.b_narrow.brc$guild_narrow <-gdf.b_narrow.brc$guild_narrow[-match(singletons, names(gdf.b_narrow.brc$guild_narrow))] ``` ```{r brcnophylo} narrowguild_model.brc <- procD.lm(Coords~mass*guild_narrow, SS.type= "II", data= gdf.b_narrow.brc, warpgrids = FALSE, iter=999, print.progress = FALSE) summary(narrowguild_model.brc) ``` ```{r brclambda} residuals.brc2<-narrowguild_model.brc$residuals residuals.pca.brc2<- prcomp(residuals.brc2) PCscores.residuals.brc2 <-residuals.pca.brc2$x[,1:13] lambdaestimate.brc2 <-transformPhylo.ML(y=PCscores.residuals.brc2, phy=gdf.b_narrow.brc$phy, model= "lambda", lowerBound = 0, upperBound = 1) lambdaestimate.brc2$Lambda ``` Lambda is not 0, so we need to correct for phylogeny in this analysis. ```{r brcphylo} transformedtree.brc2<- geiger:::rescale(gdf.b_narrow.brc$phy, "lambda", lambdaestimate.brc2$Lambda[,1]) cformvcv2 <- vcv(transformedtree.brc2) narrowguild_model.phylo<- procD.lm(Coords~mass*guild_narrow, SS.type= "II",Cov=cformvcv2, data= gdf.b_narrow.brc, warpgrids = FALSE, iter=999, print.progress = FALSE) summary(narrowguild_model.phylo) ``` Again, lambda is not 0 so we report the phylogenetically corrected version ('narrowguild_model.phylo') ## 3: Looking at the variability of allometry across families Now we will move through each of the four families with a large enough sample size. For each family, we will subset out the data and then run an ANOVA to understand how much of the variation in skull shape can be explained by foraging ecology (breeding season specific), body mass, and the interaction between the two. ```{r, message=FALSE} remove(list=ls()) load("Charadriiform_GeoMorph_Data_breeding.RData") library(Morpho) library(dplyr) library(geomorph) library(motmot) library(ggplot2) library(ape) ``` Family 1--Scolopacidae: ```{r scologdf} coords.Scolopacidae <- gdf.b$Coords[ , , gdf.b$Family == 'Scolopacidae'] mass.Scolopacidae <- gdf.b$mass[gdf.b$Family == 'Scolopacidae'] Scolopacidae.tree <- as.phylo(drop.tip(gdf.b$phy, setdiff(gdf.b$phy$tip.label, (dimnames(coords.Scolopacidae)[[3]])))) coords.Scolopacidae <- gdf.b$Coords[ , , gdf.b$Family == 'Scolopacidae'] foraging.Scolopacidae <- gdf.b$guild_narrow[gdf.b$Family == 'Scolopacidae'] cbind(Scolopacidae.tree$tip.label, names(mass.Scolopacidae), dimnames(coords.Scolopacidae)[[3]], names(foraging.Scolopacidae)) gdf.b.Scolopacidae <- geomorph.data.frame(Coords= coords.Scolopacidae, phy=Scolopacidae.tree, Weight = mass.Scolopacidae, foraging =foraging.Scolopacidae) singletons <- names(which(table(gdf.b.Scolopacidae$foraging)==1)) singletons <- names(gdf.b.Scolopacidae$foraging[match(singletons,gdf.b.Scolopacidae$foraging)]) gdf.b.Scolopacidae.trimmed<-gdf.b.Scolopacidae gdf.b.Scolopacidae.trimmed$phy <- drop.tip(gdf.b.Scolopacidae.trimmed$phy, singletons) gdf.b.Scolopacidae.trimmed$Coords<-gdf.b.Scolopacidae.trimmed$Coords[,,-match(singletons,dimnames(gdf.b.Scolopacidae.trimmed$Coords)[[3]])] gdf.b.Scolopacidae.trimmed$Weight <-gdf.b.Scolopacidae.trimmed$Weight[-match(singletons, names(gdf.b.Scolopacidae.trimmed$Weight))] gdf.b.Scolopacidae.trimmed$foraging <-gdf.b.Scolopacidae.trimmed$foraging[-match(singletons, names(gdf.b.Scolopacidae.trimmed$foraging))] ``` ```{r scoloanova} fit.Scolopacidae.trimmed <- procD.lm(Coords~Weight*foraging, SS.type= "II", data= gdf.b.Scolopacidae.trimmed, warpgrids = FALSE, iter=999, print.progress = FALSE) residuals.Scolop <-fit.Scolopacidae.trimmed$residuals residuals.Scolop <-as.matrix(residuals.Scolop) residuals.pca.Scolop<- prcomp(residuals.Scolop) summary(residuals.pca.Scolop) PCscores.residuals.Scolo <-residuals.pca.Scolop$x [,1:11] lambdaestimate.Scolop <-transformPhylo.ML(y=PCscores.residuals.Scolo, phy=gdf.b.Scolopacidae.trimmed$phy, model= "lambda") lambdaestimate.Scolop$Lambda summary(fit.Scolopacidae.trimmed) ``` Lambda here is zero, so we will report the non-phylogenetic version. But again, we can run the model with the upper estimate just to be conservative and make sure these reults don't differ much . ```{r scoloanovaupper} transformedtree.Scolop.trimmed<- geiger:::rescale(gdf.b.Scolopacidae.trimmed$phy, "lambda", lambdaestimate.Scolop$Lambda[,3]) cformvcv.Scolop <- vcv(transformedtree.Scolop.trimmed) fit.Scolopacidae.trimmed.phylo <- procD.lm(Coords~Weight*foraging , SS.type= "II",Cov= cformvcv.Scolop, data= gdf.b.Scolopacidae.trimmed, iter=999, print.progress =FALSE) summary(fit.Scolopacidae.trimmed.phylo) ``` We can also plot the small and large species morphology. ```{r, fig.show="hold", out.width="50%"} M.Scolopacidae<- mshape(gdf.b.Scolopacidae.trimmed$Coords) preds.Scolo <- shape.predictor(gdf.b.Scolopacidae.trimmed$Coords, x= (gdf.b.Scolopacidae.trimmed$Weight), Intercept = TRUE, predmin = min((gdf.b.Scolopacidae.trimmed$Weight)), predmax = max((gdf.b.Scolopacidae.trimmed$Weight))) smallspecimenmesh <-plotRefToTarget(M.Scolopacidae, preds.Scolo$predmin,method="TPS") bigspecimenmesh<- plotRefToTarget(M.Scolopacidae, preds.Scolo$predmax, method = "TPS") ``` Family 2-- Laridae (gulls, terns, and skimmers): ```{r largdf} coords.Laridae <- gdf.b$Coords[ , , gdf.b$Family == 'Laridae'] mass.Laridae <- gdf.b$mass[gdf.b$Family == 'Laridae'] Laridae.tree <- as.phylo(drop.tip(gdf.b$phy, setdiff(gdf.b$phy$tip.label, dimnames(coords.Laridae)[[3]]))) foraging.Laridae<- gdf.b$guild_narrow[gdf.b$Family == 'Laridae'] cbind(Laridae.tree$tip.label, names(mass.Laridae), dimnames(coords.Laridae )[[3]], names(foraging.Laridae)) gdf.b.Laridae<- geomorph.data.frame(Coords = coords.Laridae, phy= Laridae.tree, Weight = mass.Laridae, foraging= foraging.Laridae) singletons <- names(which(table(gdf.b.Laridae$foraging)==1)) singletons <- names(gdf.b.Laridae$foraging[match(singletons,gdf.b.Laridae$foraging)]) gdf.b.Laridae.trimmed<-gdf.b.Laridae gdf.b.Laridae.trimmed$phy <- drop.tip(gdf.b.Laridae.trimmed$phy, singletons) gdf.b.Laridae.trimmed$Coords<-gdf.b.Laridae.trimmed$Coords[,,-match(singletons,dimnames(gdf.b.Laridae.trimmed$Coords)[[3]])] gdf.b.Laridae.trimmed$Weight <-gdf.b.Laridae.trimmed$Weight[-match(singletons, names(gdf.b.Laridae.trimmed$Weight))] gdf.b.Laridae.trimmed$foraging <-gdf.b.Laridae.trimmed$foraging[-match(singletons, names(gdf.b.Laridae.trimmed$foraging))] ``` ```{r laranova} fit.Laridae.trimmed <- procD.lm(Coords~Weight*foraging, SS.type= "II", data= gdf.b.Laridae.trimmed, warpgrids = FALSE, iter=999, print.progress = FALSE) residuals.Laridae <-fit.Laridae.trimmed$residuals residuals.Laridae <-as.matrix(residuals.Laridae) residuals.pca.Laridae<- prcomp(residuals.Laridae) summary(residuals.pca.Laridae) PCscores.residuals.Laridae <-residuals.pca.Laridae$x [,1:25] lambdaestimate.Laridae <-transformPhylo.ML(y=PCscores.residuals.Laridae, phy=gdf.b.Laridae.trimmed$phy, model= "lambda") lambdaestimate.Laridae$Lambda summary(fit.Laridae.trimmed) ``` Lambda is 0 so we will report the non-phylogenetic version, but we can try the upper estimate ```{r laranovaupper} transformedtree.Laridae<- geiger:::rescale(gdf.b.Laridae.trimmed$phy, "lambda", lambdaestimate.Laridae$Lambda[,3]) cformvcv.Laridae <- vcv(transformedtree.Laridae) fit.Laridae.trimmed.phylo <- procD.lm(Coords~Weight*foraging, SS.type= "II",Cov= cformvcv.Laridae, data= gdf.b.Laridae.trimmed, warpgrids = FALSE, iter=999, print.progress = FALSE) summary(fit.Laridae.trimmed.phylo) ``` ```{r, fig.show="hold", out.width="50%"} M.Laridae <- mshape(gdf.b.Laridae.trimmed$Coords) preds.Laridae <- shape.predictor(gdf.b.Laridae.trimmed$Coords, x= (gdf.b.Laridae.trimmed$Weight), Intercept = TRUE, predmin = min((gdf.b.Laridae.trimmed$Weight)), predmax = max((gdf.b.Laridae.trimmed$Weight))) smallspecimenmesh <-plotRefToTarget(M.Laridae, preds.Laridae$predmin, method= "TPS") bigspecimenmesh<- plotRefToTarget(M.Laridae, preds.Laridae$predmax, method = "TPS") ``` Family 3-- Alcidae (Auks): ```{r alcgdf} coords.Alcids <- gdf.b$Coords[ , , gdf.b$Family == 'Alcidae'] Weights.Alcids <- gdf.b$mass[gdf.b$Family == 'Alcidae'] Alcid.tree <- as.phylo(drop.tip(gdf.b$phy, setdiff(gdf.b$phy$tip.label, dimnames(coords.Alcids)[[3]]))) gdf.b.Alcids <- geomorph.data.frame(coords = coords.Alcids, phy= Alcid.tree, Weight = Weights.Alcids) ``` Note: ALL Alcids share the same foraging mode ("Pursuit Diver, Pelagic). So we will run an ANOVA WITHOUT the foraging ecology since this has no variance. Because of this, we also do not need the type II sum of squares arguement here. Because this anlaysis does not include foraging ecology, this is an identical test to that done in our main text, but we include it here for the sake of completeness. ```{r alcidsanova} ## now linear regression for just alcids fit.Alcids <- procD.lm(coords~Weight, data= gdf.b.Alcids) #pull out the residuals, save as a matrix residuals.Alcids <-fit.Alcids$residuals residuals.Alcids <-as.matrix(residuals.Alcids) ##the residuals have too many dimensions to use as is- use PCA to reduce these residuals.pca.Alcids<- prcomp(residuals.Alcids) summary(residuals.pca.Alcids) #now save out the PC scores of the residuals- we're looking at the first 12 PCs because they together explain #about 95% of the variance PCscores.residuals.Alcids<-residuals.pca.Alcids$x [,1:11] #Now use the tree and these residuals to estimate the value of lambda lambdaestimate.Alcids <-transformPhylo.ML(y=PCscores.residuals.Alcids, phy=gdf.b.Alcids$phy, model= "lambda") lambdaestimate.Alcids$Lambda ``` Becuase lambda is NOT zero, we need to use the phylogenetically-corrected version of this ANOVA. ```{r alcidphylo} ###MLLamba is ~0.7/0.8 last decimal point may change between analyses #now rescale the tree using lambda transformedtree.Alcids<- geiger:::rescale(gdf.b.Alcids$phy, "lambda", lambdaestimate.Alcids$Lambda[,1]) ###pull out the variance covariance matrix from the rescaled tree and re-run the ANOVA using this covariance matrix cformvcv.Alcids <- vcv(transformedtree.Alcids) fit.Alcids_phylo <- procD.lm(coords~Weight,Cov= cformvcv.Alcids, data= gdf.b.Alcids, iter=999, print.progress = FALSE) summary(fit.Alcids_phylo) ``` ```{r, fig.show="hold", out.width="50%"} M.Alcids <- mshape(coords.Alcids) preds.Alcids <- shape.predictor(gdf.b.Alcids$coords, x= gdf.b.Alcids$Weight, Intercept = TRUE, predmin = min((gdf.b.Alcids$Weight)), predmax = max((gdf.b.Alcids$Weight))) smallspecimenmesh <-plotRefToTarget(M.Alcids, preds.Alcids$predmin, method= "TPS") bigspecimenmesh<- plotRefToTarget(M.Alcids, preds.Alcids$predmax, method = "TPS") ``` Family 4 -- Charadriidae (plovers, lapwings, and relatives): ```{r chargdf} coords.Chardriidae <- gdf.b$Coords[ , , gdf.b$Family == 'Charadriidae'] Weight.Charadriidae <- gdf.b$mass[gdf.b$Family == 'Charadriidae'] foraging.Charadriidae <- gdf.b$guild_narrow[gdf.b$Family == 'Charadriidae'] Charadriidae.tree <- as.phylo(drop.tip(gdf.b$phy, setdiff(gdf.b$phy$tip.label, dimnames(coords.Chardriidae)[[3]]))) gdf.b.Charadriidae<- geomorph.data.frame(coords=coords.Chardriidae, phy= Charadriidae.tree, Weight = Weight.Charadriidae, foraging= foraging.Charadriidae) singletons <- names(which(table(gdf.b.Charadriidae$foraging)==1)) singletons <- names(gdf.b.Charadriidae$foraging[match(singletons,gdf.b.Charadriidae$foraging)]) gdf.b.Charadriidae.trimmed<-gdf.b.Charadriidae gdf.b.Charadriidae.trimmed$phy <- drop.tip(gdf.b.Charadriidae.trimmed$phy, singletons) gdf.b.Charadriidae.trimmed$coords<-gdf.b.Charadriidae.trimmed$coords[,,-match(singletons,dimnames(gdf.b.Charadriidae.trimmed$coords)[[3]])] gdf.b.Charadriidae.trimmed$Weight <-gdf.b.Charadriidae.trimmed$Weight[-match(singletons, names(gdf.b.Charadriidae.trimmed$Weight))] gdf.b.Charadriidae.trimmed$foraging <-gdf.b.Charadriidae.trimmed$foraging[-match(singletons, names(gdf.b.Charadriidae.trimmed$foraging))] ``` ```{r charanova} fit.Char.trimmed <- procD.lm(coords~Weight*foraging, SS.type= "II", data= gdf.b.Charadriidae.trimmed, warpgrids = FALSE, iter=999, print.progress = FALSE) residuals.Char <-fit.Char.trimmed$residuals residuals.Char <-as.matrix(residuals.Char) residuals.pca.Char<- prcomp(residuals.Char) summary(residuals.pca.Char) PCscores.residuals.Char<-residuals.pca.Char$x [,1:17] lambdaestimate.Char <-transformPhylo.ML(y=PCscores.residuals.Char, phy=gdf.b.Charadriidae.trimmed$phy, model= "lambda") lambdaestimate.Char$Lambda summary(fit.Char.trimmed) ``` Lambda is zero, so we will report the non-phylogenetically corrected version. But, again, we can run this model using the upper estimate of lambda just to be safe. ```{r charupper} transformedtree.Char<- geiger:::rescale(gdf.b.Charadriidae.trimmed$phy, "lambda", lambdaestimate.Char$Lambda[,3]) cformvcv.Char<- vcv(transformedtree.Char) fit.Char_phylo <- procD.lm(coords~Weight*foraging, SS.type= "II",Cov= cformvcv.Char, data= gdf.b.Charadriidae.trimmed, warpgrids = FALSE, iter=999, print.progress = FALSE) summary(fit.Char_phylo) ``` ```{r, fig.show="hold", out.width="50%"} M.Char <- mshape(gdf.b.Charadriidae.trimmed$coords) preds.Char <- shape.predictor(gdf.b.Charadriidae.trimmed$coords, x= gdf.b.Charadriidae.trimmed$Weight, Intercept = TRUE, predmin = min((gdf.b.Charadriidae.trimmed$Weight)), predmax = max((gdf.b.Charadriidae.trimmed$Weight))) smallspecimenmesh <-plotRefToTarget(M.Char, preds.Char$predmin, method= "TPS") bigspecimenmesh<- plotRefToTarget(M.Char, preds.Char$predmax, method = "TPS") ```