Package 'geomorph'

Description

Functions in this package allow one to read, manipulate, and digitize landmark data; generate shape variables via Procrustes analysis for points, curves and surface data, perform statistical analyses of shape variation and covariation, and provide graphical depictions of shapes and patterns of shape variation.

geomorph TOC

geomorph-package

Author(s)

Dean C. Adams, Michael Collyer, Antigoni Kaliontzopoulou, and Erica Baken

Description

Convert a matrix of landmark coordinates into a three-dimensional array

Usage

arrayspecs(A, p, k, sep = NULL)

Arguments

Details

This function converts a matrix of landmark coordinates into a 3D array (p x k x n), which is the required input format for many functions in geomorph. The input matrix can be arranged such that the coordinates of each landmark are found on a separate row, or that each row contains all landmark coordinates for a single specimen.

Value

Function returns a 3D array (p x k x n), where p is the number of landmark points, k is the number of landmark dimensions (2 or 3), and n is the number of specimens. The third dimension of this array contains names for each specimen if specified in the original input matrix.

Author(s)

Dean Adams & Mike Collyer

See Also

two.d.array

Examples

## Not run: 

x <- matrix(rnorm(18), nrow = 3)  # Random triangles (all coordinates on same 
 # row for each triangle)
arrayspecs(x, 3, 2) 
 
x2 <- matrix(rnorm(18), ncol = 2) # Random triangles (each landmark on its 
# own row)
arrayspecs(x2, 3, 2)

## End(Not run)

Description

Function performs an analysis of directional and fluctuating asymmetry for bilaterally symmetric objects

Usage

bilat.symmetry(
  A,
  ind = NULL,
  side = NULL,
  replicate = NULL,
  object.sym = FALSE,
  land.pairs = NULL,
  data = NULL,
  iter = 999,
  seed = NULL,
  RRPP = TRUE,
  SS.type = c("I", "II", "III"),
  turbo = TRUE,
  Parallel = FALSE,
  print.progress = TRUE,
  ...
)

Arguments

Details

The function quantifies components of shape variation for a set of specimens as described by their patterns of symmetry and asymmetry. Here, shape variation is decomposed into variation among individuals, variation among sides (directional asymmetry), and variation due to an individual x side interaction (fluctuating symmetry). These components are then statistically evaluated using Procrustes ANOVA. Statistical assessment of model effects for shape variation is accomplished using permutation procedures. Methods for both matching symmetry and object symmetry can be implemented. Matching symmetry is when each object contains mirrored pairs of structures (e.g., right and left hands) while object symmetry is when a single object is symmetric about a midline (e.g., right and left sides of human faces). Details on general approaches for the study of symmetry in geometric morphometrics may be found in: Mardia et al. 2000; Klingenberg et al. 2002.

As input, the function receives either A 3D array (p x k x n) containing raw landmarks (requiring GPA to be performed) or a gpagen object (if GPA has been previously performed) or a geomorphShapes object. If one wishes to incorporate semilandmarks, GPA can either be performed first using gpagen, or within bilat.symmetry by passing adequate GPA arguments (i.e. curves, surfaces, ProcD etc, see gpagen. If a geomorphShapes object is provided, semilandmarks are automatically identified and slid during GPA. For object.sym = FALSE, landmarks should be of dimension (p x k x 2n), as each specimen is represented by both left and right configurations.

Analyses of symmetry for matched pairs of objects is implemented when object.sym = FALSE. Here, a 3D array [p x k x 2n] contains the landmark coordinates for all pairs of structures (2 structures for each of n specimens). Because the two sets of structures are on opposite sides, they represent mirror images, and one set must be reflected prior to the analysis to allow landmark correspondence. IT IS ASSUMED THAT THE USER HAS DONE THIS PRIOR TO PERFORMING THE SYMMETRY ANALYSIS. Reflecting a set of specimens may be accomplished by multiplying one coordinate dimension by '-1' for these structures (either the x-, the y-, or the z-dimension). A vector containing information on individuals and sides must also be supplied. Replicates of each specimen may also be included in the dataset, and when specified will be used as measurement error (see Klingenberg and McIntyre 1998).

Analyses of object symmetry is implemented when object.sym = TRUE. Here, a 3D array [p x k x n] contains the landmark coordinates for all n specimens. To obtain information about asymmetry, the function generates a second set of objects by reflecting them about one of their coordinate axes. The landmarks across the line of symmetry are then relabeled to obtain landmark correspondence. The user must supply a list of landmark pairs. A vector containing information on individuals must also be supplied. Replicates of each specimen may also be included in the dataset, and when specified will be used as measurement error.

The function also provides individual measures of unsigned asymmetry, calculated as the Procrustes distance between the right and left element (for paired structures, as detailed in Klingenberg and McIntyre 1998) or side of the structure (for object symmetry, following Lazić et al 2015). The computational difference between the two approaches consists in that, for object symmetry, only paired landmarks are considered, excluding the landmarks of the midline.

Notes for geomorph 3.0

Compared to older versions of geomorph, some results can be expected to be slightly different. Starting with geomorph 3.0, results use only type I sums of squares (SS) with either full randomization of raw shape values or RRPP (preferred with nested terms) for analysis of variance (ANOVA). Older versions used a combination of parametric and non-parametric results, as well as a combination of type I and type III SS. While analytical conclusions should be consistent (i.e., "significance" of effects is the same), these updates maintain consistency in analytical philosophy. This change will require longer computation time for large datasets, but the trade-off allows users to have more flexibility and eliminates combining disparate analytical philosophies.

Note also that significance of terms in the model are found by comparing F-values for each term to those obtained via permutation. F-ratios and df are not strictly necessary (a ratio of SS would suffice), but they are reported as is standard for anova tables. Additionally, users will notice that the df reported are based on the number of observations rather than a combination of objects * coordinates * dimensions, as is sometimes found in morphometric studies of symmetry. However, this change has no effect on hypothesis testing, as only SS vary among permutations (df, coordinates, and dimensions are constants).

The generic functions, print, summary, and plot all work with bilat.symmetry.

Value

An object of class "bilat.symmetry" returns a list of the following

Author(s)

Dean Adams, Michael Collyer and Antigoni Kaliontzopoulou

References

Klingenberg, C.P. and G.S. McIntyre. 1998. Quantitative genetics of geometric shape in the mouse mandible. Evolution. 55:2342-2352.

Mardia, K.V., F.L. Bookstein, and I.J. Moreton. 2000. Statistical assessment of bilateral symmetry of shapes. Biometrika. 87:285-300.

Klingenberg, C.P., M. Barluenga, and A. Meyer. 2002. Shape analysis of symmetric structures: quantifying variation among individuals and asymmetry. Evolution. 56:1909-1920.

Lazić, M. M., M. A. Carretero, J. Crnobrnja-Isailović, and A. Kaliontzopoulou. 2015. Effects of environmental disturbance on phenotypic variation: an integrated assessment of canalization, developmental stability, modularity, and allometry in lizard head shape. The American Naturalist 185:44–58.

Examples

## Not run: 

#Example of matching symmetry
data(mosquito)
gdf <- geomorph.data.frame(wingshape = mosquito$wingshape, 
ind = mosquito$ind, 
side = mosquito$side,
replicate = mosquito$replicate)
mosquito.sym <- bilat.symmetry(A = wingshape, ind = ind, side = side,
replicate = replicate, object.sym = FALSE, RRPP = TRUE, 
data = gdf)
summary(mosquito.sym)
plot(mosquito.sym, warpgrids = TRUE)
mosquito.sym$shape.anova # extract just the anova table on shape

# Previous example, performing GPA first
Y.gpa <- gpagen(mosquito$wingshape)
mosquito.sym2 <- bilat.symmetry(A = Y.gpa, ind = ind, side = side,
replicate = replicate, object.sym = FALSE, RRPP = TRUE, 
data = gdf)
summary(mosquito.sym2)
summary(mosquito.sym) # same results

#Example of object symmetry

data(lizards)
gdf <- geomorph.data.frame(shape = lizards$coords, 
ind = lizards$ind, 
replicate = lizards$rep)
liz.sym <- bilat.symmetry(A = shape, ind = ind, rep = rep, 
object.sym = TRUE, 
land.pairs = lizards$lm.pairs, data = gdf, RRPP = TRUE)
summary(liz.sym)

# Example of object symmetry in 3D and including semilandmarks

data(scallops)
gdf <- geomorph.data.frame(shape = scallops$coorddata, 
ind = scallops$ind)
scallop.sym <- bilat.symmetry(A = shape, ind = ind, 
object.sym = TRUE, 
curves= scallops$curvslide, surfaces = scallops$surfslide,
land.pairs=scallops$land.pairs, data = gdf, RRPP = TRUE)
summary(scallop.sym)
# NOTE one can also: plot(scallop.sym, warpgrids = TRUE, mesh = NULL)
# NOTE one can also: scallop.sym$data.type # recall the symmetry type

## End(Not run)

Description

An interactive function to build a template for the digitization across specimens of three dimensional (3D) surface sliding semilandmarks. Input for the function is either a matrix of vertex coordinates defining a 3D surface object or a mesh3d object as obtained from read.ply.

Usage

buildtemplate(spec, fixed, center = TRUE, surface.sliders, ptsize = 1)

Arguments

Details

Function constructs a template of surface slider semilandmarks. If desired, the user can simultaneously digitize the fixed points (see digitizing below), however these may have been previously digitized separately using digit.fixed.

The function finds surface semilandmarks, chosen automatically from the mesh at roughly equal distances, using the nearest-neighbor algorithm outlined in Gunz et al. (2005) and Mitteroecker and Gunz (2009).

The set of fixed and surface slider landmarks are saved in the working directory as a txt file named "template". This file will then be used to extract a set of similarly numbered surface semilandmarks on subsequent specimens using the function digitsurface. Because template matching is based on the correspondence of fixed landmark points in the template and the target specimen, a minimum of four fixed landmarks must be used. However, to ensure a strong match between the scan and the template, it is recommended that a higher number of fixed points is used.

For more details see the vignette: vignette("geomorph.digitize3D").

NOTE: Function centers the mesh before digitizing by default (center = TRUE). If one chooses not to center, the specimen may be difficult to manipulate in the rgl window.

Digitizing

Digitizing of fixed landmarks is interactive. Once a point is selected, the user is asked if the system should keep or discard the selection (y/n). If "y", the user is asked to continue to select the next landmark. If "n" the removes the last chosen landmark, and the user is asked to select it again. This can be repeated until the user is comfortable with the landmark chosen.

To digitize with a standard 3-button (PC):

1. the RIGHT mouse button (primary) to select points to be digitized,

2. the LEFT mouse button (secondary) is used to rotate mesh,

3. the mouse SCROLLER (third/middle) is used to zoom in and out.

NOTE: Digitizing functions on MACINTOSH computers using a standard 3-button mice works as specified. Macs using platform specific single button mice, XQuartz must be configured: go to Preferences > Input > tick "Emulate three button mouse":

1. press button to rotate 3D mesh,

2. press button while pressing COMMAND key to select vertex to be used as a landmark),

3. press button while pressing OPTION key to adjust mesh perspective.

4. the mouse SCROLLER or trackpad two finger scroll is used to zoom in an out.

NOTE: there is no pan (translate) functionality in rgl library for all platforms at this time. The template can be edited using function editTemplate.

AUTO mode

The function as described above (for interactive mode) calls digit.fixed, prompting the user to select fixed landmarks in the rgl window. However if the user has digitized these fixed landmark elsewhere (e.g., in other software), then the input for parameter 'fixed' can be a p-x-k matrix of 3D coordinates. In this case, the function will automatically use these landmarks to build the template of sliding semilandmarks.

Value

The function writes to the working directory three files: an NTS file with the name of the specimen and .nts suffix containing the landmark coordinates, "template.txt" containing the same coordinates for use with the function digitsurface, and "surfslide.csv", a file containing the address of the landmarks defined as "surface sliders" for use with gpagen. The function also returns to the console an n x 3 matrix containing the x,y,z coordinates of the digitized landmarks.

Author(s)

Erik Otarola-Castillo & Emma Sherratt

References

Gunz P, Mitteroecker P, & Bookstein FJ (2005) Semilandmarks in Three Dimensions. Modern Morphometrics in Physical Anthropology, ed Slice DE (Springer-Verlag, New York), pp 73-98.

Mitteroecker P & Gunz P (2009) Advances in Geometric Morphometrics. Evolutionary Biology 36(2):235-247.

See Also

read.ply

digit.fixed

digitsurface

Description

Combine separate landmark configurations (subsets) into one landmark set

Usage

combine.subsets(
  ...,
  gpa = TRUE,
  CS.sets = NULL,
  norm.CS = FALSE,
  weights = NULL
)

Arguments

Details

This function combines landmark configurations (either landmarks requiring GPA or Procrustes shape variables following GPA) to create a different morphological data set. This might be of interest, for example, if one has landmarks digitized on separate images collected from the same organisms. (In the examples below, configurations for heads and tails of larval salamanders were collected separately from images taken on the same individuals.) An attempt is made to scale configurations by their relative centroid sizes, following the procedure in Davis et al. (2016); i.e., landmark coordinates are multiplied by CSi/sqrt(CSi^2 + CSj^2 + ...) before combining them, so that resulting combinations of landmarks are scaled to unit centroid size. This is only possible if GPA is performed on landmarks (gpa = TRUE) or centroid sizes are provided as an argument. Objects of class gpagen can be used rather than original landmarks (recommended, especially if curves or surface sliding semilandmarks are used, as different arguments cannot be passed onto onto separate GPAs via this function).

The procedure of Davis et al. (2016) is an extension of the "separate subsets" method of Adams (1999) for articulated structures.

Value

An object of class combined.set is a list containing the following

Author(s)

Michael Collyer

References

Davis, M.A., M.R. Douglas, M.L. Collyer, & M.E. Douglas, M. E. 2016. Deconstructing a species-complex: geometric morphometric and molecular analyses define species in the Western Rattlesnake (Crotalus viridis). PLoS one, 11(1), e0146166.

Adams, D.C. 1999. Methods for shape analysis of landmark data from articulated structures. Evolutionary Ecology Research. 1:959-970.

Dryden, I.L. and K.V Mardia. 2016. Statistical shape analysis, with applications in R: Second edition.

Collyer, M.L., M.A. Davis, and D.C. Adams. 2020. Making heads or tails of combined landmark configurations in geometric morphometric data. Evolutionary Biology. 47:193-205.

Examples

## Not run: 

data(larvalMorph) 
 head.gpa <- gpagen(larvalMorph$headcoords, 
   curves = larvalMorph$head.sliders)
 tail.gpa <- gpagen(larvalMorph$tailcoords, 
 curves = larvalMorph$tail.sliders)

# Combine original data without GPA (plot to see relative size of  
# heads and tails)

 all.lm <- combine.subsets(head = larvalMorph$headcoords,
 tail = larvalMorph$tailcoords, gpa = FALSE, CS.sets = NULL)
 plotAllSpecimens((all.lm$coords))
 
 # Combine with GPA and relative centroid size
 
comb.lm <- combine.subsets(head = head.gpa, tail = tail.gpa, gpa = TRUE)
summary(comb.lm)

(configurations are actual relative size)
comb.lm$coords[,,1]

# Plot all specimens and just first specimen and color code landmarks 
par(mfrow = c(1,2))
plotAllSpecimens(comb.lm$coords)
plot(comb.lm$coords[,,1], pch = 21, bg = c(rep(1,26), 
rep(2,64)), asp = 1)

# Override relative centroid size

comb.lm <- combine.subsets(head = head.gpa$coords, 
tail = tail.gpa$coords, gpa = FALSE, CS.sets = NULL)
par(mfrow = c(1,2))
plotAllSpecimens(comb.lm$coords)
plot(comb.lm$coords[,,1], pch = 21, bg = c(rep(1,26), 
rep(2,64)), asp = 1)

# Note the head is as large as the tail, which is quite unnatural.

## Normalizing centroid size

comb.lm <- combine.subsets(head = head.gpa, 
tail = tail.gpa, gpa = TRUE, norm.CS = TRUE)
summary(comb.lm)
par(mfrow = c(1,2))
plotAllSpecimens(comb.lm$coords)
plot(comb.lm$coords[,,1], pch = 21, bg = c(rep(1,26), 
rep(2,64)), asp = 1)
par(mfrow = c(1,1))

# Note that the head is too large, compared to a real specimen.  
# This option focuses on average distance of points to centroid, 
# but ignores the number of landmarks.  
# Consequently,the density of landmarks in the head and tail are 
# irrelevant and the head size is inflated because of the fewer 
# landmarks in the configuration.

## Weighting centroid size

comb.lm <- combine.subsets(head = head.gpa, 
tail = tail.gpa, gpa = TRUE, norm.CS = FALSE, weights = c(0.3, 0.7))
summary(comb.lm)
par(mfrow = c(1,2))
plotAllSpecimens(comb.lm$coords)
plot(comb.lm$coords[,,1], pch = 21, bg = c(rep(1,26), 
rep(2,64)), asp = 1)
par(mfrow = c(1,1))

# Note that the head is way too small, compared to a real specimen.  
# This option allows one to dictate the relative sizes of subsets
#  as portions of the combined set.  An option like this should be 
# used with caution, but can help overcome issues caused by landmark 
# density.

## End(Not run)

Description

Function performs an analysis to compare the effect sizes of two or more CR effects

Usage

compare.CR(..., CR.null = TRUE, two.tailed = TRUE)

Arguments

Details

The function statistically compares the effect sizes of two or more CR analyses. Typically, this function might be used to compare levels of modularity between two or more samples, each measuring the degree of morphological modularity in each. Alternatively, the approach can compare the degree of modular signal as expressed by alternative modular hypotheses for the same dataset.

The analysis calculates effect sizes as standard deviates, z, and performs two-sample z-tests, using the pooled standard error from the sampling distributions of the CR analyses. The method follows that of Adams and Collyer (2019) used to compare patterns of modularity across datasets.

To use this function, simply perform modularity.test, or phylo.modularity on as many samples or alternative modular hypotheses as desired. Any number of objects of class CR can be input. For the case of the latter, one may wish to include the null hypothesis of no modularity (i.e., that all variables belong to a single module). For this, the CR.null = TRUE option should be specified. Finally, one may perform the comparison as either a one-tailed or a two-tailed (default) test.

Value

An object of class compare.CR, returns a list of the following

Author(s)

Dean Adams and Michael Collyer

References

Adams, D.C. and M.L. Collyer. 2019. Comparing the strength of modular signal, and evaluating alternative modular hypotheses, using covariance ratio effect sizes for morphometric data. Evolution. 73:2352-2367.

Examples

## Not run: 

# Example 1: Compare modular signal across datasets
 
 data(pupfish) 
 Y.gpa<-gpagen(pupfish$coords, print.progress = FALSE)    #GPA-alignment   
 
## landmarks on the body and operculum
 land.gps<-rep('a',56); land.gps[39:48]<-'b'

 group <- factor(paste(pupfish$Pop, pupfish$Sex, sep = "."))
 levels(group)

 coords.gp <- coords.subset(Y.gpa$coords, group)

 modul.tests <- Map(function(x) modularity.test(x, land.gps,iter=999, 
 print.progress = FALSE), coords.gp) 
         
# the map function performs the integration test on each 3D array 
# in the lists provided

  modul.tests$Marsh.F
  modul.tests$Marsh.M
  modul.tests$Sinkhole.F
  modul.tests$Sinkhole.M

 group.Z <- compare.CR(modul.tests, CR.null = FALSE)
 summary(group.Z)

# Example 2: Compare alternative modular hypotheses

# 3 module hypothesis (tail now a module)
 land.gps3 <- rep('a',56); land.gps3[39:48]<-'b'; 
 land.gps3[c(6:9,28:38)] <- 'c' 
   
# 4 module hypothesis (eye now a module)
 land.gps4 <- rep('a',56); land.gps4[39:48]<-'b'; 
 land.gps4[c(6:9,28:38)] <- 'c'; 
  land.gps4[c(10,49:56)] <- 'd'  

 m3.test <- modularity.test(coords.gp$Marsh.F,land.gps3, 
 print.progress = FALSE)
 m4.test <- modularity.test(coords.gp$Marsh.F,land.gps4, 
 print.progress = FALSE)

 model.Z <- compare.CR(modul.tests$Marsh.F,m3.test,m4.test, 
 CR.null = TRUE)
 summary(model.Z)

## End(Not run)

Description

Function calculates net rates of shape evolution for two or more groups of species on a phylogeny from a set of Procrustes-aligned specimens

Usage

compare.evol.rates(
  A,
  phy,
  gp,
  iter = 999,
  seed = NULL,
  method = c("permutation", "simulation"),
  print.progress = TRUE
)

Arguments

Details

The function compares net rates of morphological evolution for two or more groups of species on a phylogeny, under a Brownian motion model of evolution. It is assumed that the landmarks have previously been aligned using Generalized Procrustes Analysis (GPA) [e.g., with gpagen]. The approach is based on the outer-product matrix of between species differences in morphospace after phylogenetic transformation (Adams 2014). From the data the net rate of shape evolution for each group in the multi-dimensional space is calculated, and a ratio of rates is obtained. If three or more groups of species are used, the ratio of the maximum to minimum rate is used as a test statistic (see Adams 2014). The function can be used with univariate data (i.e. centroid size) if imported as matrix with rownames giving the taxa names.

The generic functions, print, summary, and plot all work with compare.evol.rates. The generic function, plot, produces a histogram of random rate-ratios associated with the resampling procedure.

Notes for geomorph 3.0.4 and subsequent versions

Significance testing is now accomplished in one of two ways. First, phylogenetic simulation may be used, in which tips data are obtained under Brownian motion using a common evolutionary rate pattern for all species on the phylogeny. Specifically, the common evolutionary rate matrix for all species is used, with the multi-dimensional rate used along the diagonal elements (see Denton and Adams 2015). This procedure is more general than the original simulation procedure, and retains the desirable statistical properties of earlier methods, and under a wider array of data types. Second, significance may be accomplished via permutation, where data values at the tips are permuted relative to the (see Adams and Collyer 2018). This procedure is shown to retain all appropriate statistical properties, including rotation-invariance of significance levels (see results of Adams and Collyer 2018). In addition, a multivariate effect size describing the strength of the effect is estimated from the empirically-generated sampling distribution (see details in Adams and Collyer 2019). Values from these distributions are log-transformed prior to effect size estimation, to assure normally distributed data.

Value

An object of class "evolrate" returns a list with the following components:

Author(s)

Dean Adams & Emma Sherratt

References

Adams, D.C. 2014. Quantifying and comparing phylogenetic evolutionary rates for shape and other high-dimensional phenotypic data. Syst. Biol. 63:166-177.

Denton, J.S.S., and D.C. Adams. 2015. A new phylogenetic test for comparing multiple high-dimensional evolutionary rates suggests interplay of evolutionary rates and modularity in lanternfishes (Myctophiformes; Myctophidae). Evolution. 69:2425-2440.

Adams, D.C. and M.L. Collyer. 2018. Multivariate comparative methods: evaluations, comparisons, and recommendations. Systematic Biology. 67:14-31.

Adams, D.C. and M.L. Collyer. 2019. Comparing the strength of modular signal, and evaluating alternative modular hypotheses, using covariance ratio effect sizes with morphometric data. Evolution. 73:2352-2367.

Examples

## Not run: 

data(plethspecies) 
Y.gpa <- gpagen(plethspecies$land)    #GPA-alignment    
 gp.end <- factor(c(0,0,1,0,0,1,1,0,0))  #endangered species vs. rest
 names(gp.end) <- plethspecies$phy$tip

ER<-compare.evol.rates(A = Y.gpa$coords, phy = plethspecies$phy,
  method = "simulation", gp = gp.end)
summary(ER)
plot(ER)

## End(Not run)

Description

Function calculates net rates of shape evolution for two or more multi-dimensional traits on a phylogeny from a set of Procrustes shape variables

Usage

compare.multi.evol.rates(
  A,
  gp,
  phy,
  Subset = TRUE,
  iter = 999,
  seed = NULL,
  print.progress = TRUE
)

Arguments

Details

The function compares net rates of morphological evolution for two or more multi-dimensional traits on a phylogeny, under a Brownian motion model of evolution following the procedure of Denton and Adams (2015). It is assumed that the landmarks for all traits have previously been aligned using Generalized Procrustes Analysis (GPA) [e.g., with gpagen]. The approach calculates multivariate net evolutionary rates found from the outer-product matrix of between species differences in morphospace after phylogenetic transformation (sensu Adams 2014). From the data the net rate of shape evolution for each multi-dimensional trait is calculated, and a ratio of rates is obtained. If three or more traits are used, the ratio of the maximum to minimum rate is used as a test statistic (see Denton and Adams 2015). Significance testing is accomplished by phylogenetic simulation in which tips data are obtained under Brownian motion using a an evolutionary rate matrix for all traits, which contains a common rate for all trait dimensions (Denton and Adams 2015). If three or more traits are used, pairwise p-values are also returned. In addition, a multivariate effect size describing the strength of the effect is estimated from the empirically-generated sampling distribution (see details in Adams and Collyer 2019). Values from these distributions are log-transformed prior to effect size estimation, to assure normally distributed data.

The shape data may be input as either a 3D array (p x k x n) containing Procrustes shape variables for a set of species, or as a matrix (n x [p x k]) whose rows correspond to each species. In both cases, species names must be provided as rownames (for a matrix) or as the names of the third dimension of the array. Landmark groups for each trait are then specified by a factor array designating which landmark belongs to which trait. Additionally, if the method is to be used with other data (i.e., a set of length measurements), the input A should be a matrix of n rows of species and p columns of variables. In this case, the grouping factor should have each variable assigned to a trait group.

Comparisons of net evolutionary rates between traits may be accomplished in one of two ways. First, if the traits are are part of a single shape that was subjected to a single Procrustes superimposition (i.e., they are subsets of landmarks in the configuration), then the procedure is performed without alteration as described above. However, if the shapes are derived from different structures (shapes) that were superimposed separately, then the estimates of the rates must take the difference in the number of trait dimensions into account (see discussion in Denton and Adams 2015). This option is identified by selecting Subset = FALSE.

With compare.multi.evol.rates, the generic functions print, summary, and plot all work. The generic function, plot, produces a histogram of random rate-ratios associated with the resampling procedure.

Value

An object of class "evolrate" returns a list of the following:

Author(s)

Dean Adams (used in some internal computations)

References

Adams, D.C. 2014. Quantifying and comparing phylogenetic evolutionary rates for shape and other high-dimensional phenotypic data. Syst. Biol. 63:166-177.

Denton, J.S.S., and D.C. Adams. 2015. A new phylogenetic test for comparing multiple high-dimensional evolutionary rates suggests interplay of evolutionary rates and modularity in lanternfishes (Myctophiformes; Myctophidae). Evolution. 69:2425-2440.

Adams, D.C. and M.L. Collyer. 2019. Comparing the strength of modular signal, and evaluating alternative modular hypotheses, using covariance ratio effect sizes with morphometric data. Evolution. 73:2352-2367.

Examples

## Not run: 

data(plethspecies) 
Y.gpa <- gpagen(plethspecies$land)    #GPA-alignment    
land.gp <- c("A","A","A","A","A","B","B","B","B","B","B")  
    #mandible and cranium subsets

EMR <- compare.multi.evol.rates(A = Y.gpa$coords, gp = land.gp, 
    Subset = TRUE, phy = plethspecies$phy)
summary(EMR)
plot(EMR)

## End(Not run)

Description

Function performs an analysis to compare the effect sizes of two or more phylogenetic effect sizes

Usage

compare.physignal.z(..., two.tailed = TRUE)

Arguments

Details

The function statistically compares the effect sizes of two or more physignal.z analyses. This can be performed on different traits from the same tree, same or different traits from different trees, or modules of landmark configurations.

To use this function, perform physignal.z on as many samples as desired. Any number of objects of class physignal.z can be input. Note that some values of Z can be NaN, if the scaling parameter, lambda, is optimized at 0. In these cases, the standard error is also 0, and pairwise comparisons might not make sense.

Value

An object of class compare.physignal.z, returns a list of the following

Author(s)

Michael Collyer

References

Collyer, M.L., E.K. Baken, & D.C. Adams. 2022. A standardized effect size for evaluating and comparing the strength of phylogenetic signal. Methods in Ecology and Evolution. 13:367-382.

Examples

## Not run: 

# Example: Compare phylogenetic signal of head components in Plethodon

data(plethspecies) 
Y.gpa <- gpagen(plethspecies$land)    #GPA-alignment

## landmarks of the jaw and cranium
jaw <- 1:5
cranium <- 6:11

PS.jaw <- physignal.z(A = Y.gpa$coords[jaw,,], phy = plethspecies$phy, 
lambda = "front", PAC.no = 7, iter=999)

PS.cranium <- physignal.z(A = Y.gpa$coords[cranium,,], phy = plethspecies$phy, 
lambda = "front", PAC.no = 7, iter=999)

PS.list <-list(PS.jaw, PS.cranium)
names(PS.list) <- c("jaw", "cranium")

PS.Z <- compare.physignal.z(PS.list)
summary(PS.Z)

## End(Not run)

Description

Function performs an analysis to compare the effect sizes of two or more PLS effects

Usage

compare.pls(..., two.tailed = TRUE)

Arguments

Details

The function statistically compares the effect sizes of two or more PLS analyses. Typically, this function might be used to compare levels of integration between two or more samples, each measuring morphological integration between different modules. In such cases, the PLS correlation coefficient, r, is not a good measure of integration effect, as its expected value is dependent on both the number of specimens and number of variables (Adams and Collyer 2016). This analysis calculates effect sizes as standard deviates, z, and performs two-sample z-tests, using the pooled standard error from the sampling distributions of the PLS analyses.

To use this function, perform two.b.pls, integration.test, or phylo.integration on as many samples as desired. Any number of objects of class pls can be input.

Similar versions of this function will be designed for alternative test statistics, in the future.

Notes for geomorph 3.0.4 and subsequent versions

Compared to previous versions of geomorph, users might notice differences in effect sizes. Previous versions used z-scores calculated with expected values of statistics from null hypotheses (sensu Collyer et al. 2015); however Adams and Collyer (2016) showed that expected values for some statistics can vary with sample size and variable number, and recommended finding the expected value, empirically, as the mean from the set of random outcomes. Geomorph 3.0.4 and subsequent versions now center z-scores on their empirically estimated expected values and where appropriate, log-transform values to assure statistics are normally distributed. This can result in negative effect sizes, when statistics are smaller than expected compared to the average random outcome. For ANOVA-based functions, the option to choose among different statistics to measure effect size is now a function argument.

Value

An object of class compare.pls, returns a list of the following

Author(s)

Michael Collyer

References

Collyer, M.L., D.J. Sekora, and D.C. Adams. 2015. A method for analysis of phenotypic change for phenotypes described by high-dimensional data. Heredity. 115:357-365.

Adams, D.C. and M.L. Collyer. 2016. On the comparison of the strength of morphological integration across morphometric datasets. Evolution. 70:2623-2631.

Examples

## Not run: 
# Example of comparative morphological integration between pupfish head 
# and body shapes

 data(pupfish) # GPA previously performed
  
 group <- factor(paste(pupfish$Pop, pupfish$Sex, sep = "."))
 levels(group)
  
 tail.LM <- c(1:3, 5:9, 18:38)
 head.LM <- (1:56)[-tail.LM]

 tail.coords <- pupfish$coords[tail.LM,,]
 head.coords <- pupfish$coords[head.LM,,]
 
 # Subset 3D array by group, returning a list of 3D arrays
 tail.coords.gp <- coords.subset(tail.coords, group)
 head.coords.gp <- coords.subset(head.coords, group)

 integ.tests <- Map(function(x,y) integration.test(x, y, iter=499, 
 print.progress = FALSE), head.coords.gp, tail.coords.gp)
# the map function performs the integration test on each 3D array in 
# the lists provided

 integ.tests$Marsh.F
 integ.tests$Marsh.M
 integ.tests$Sinkhole.F
 integ.tests$Sinkhole.M

 group.Z <- compare.pls(integ.tests)
 summary(group.Z)

 # Sexual dimorphism in morphological integration in one population
 # but not the other

 # can also list different PLS analyses, separately

 compare.pls(MF = integ.tests$Marsh.F, MM = integ.tests$Marsh.M)

## End(Not run)

Description

Function performs an analysis to compare the effect sizes of two or more Vrel effects

Usage

compare.ZVrel(..., two.tailed = TRUE)

Arguments

Details

The function statistically compares the effect sizes of two or more Vrel analyses. Typically, this function might be used to compare the strength of integration in one dataset as compared with another (see Conaway and Adams 2022).

The analysis performs two-sample z-tests based on effect sizes (Z-scores) of Vrel. The method follows that of Conaway Adams (2022) used to compare the strength of integration across datasets.

To use this function, simply perform integration.Vrel on as many samples or as desired. Any number of objects of class rel.eig can be input. One may perform the comparison as either a one-tailed or a two-tailed (default) test.

Value

An object of class compare.rel.eig, returns a list of the following

Author(s)

Dean Adams

References

Conaway, M.A., and D.C. Adams. 2022. An effect size for comparing the strength of morphological integration across studies. Evolution. 76: 2244-2259.

Examples

## Not run: 
 data("plethodon")
 Y.gpa <- gpagen(plethodon$land)
 
 coords.gp <- coords.subset(Y.gpa$coords, plethodon$species)
 Vrel.gp <- Map(function(x) integration.Vrel(x), coords.gp) 
 
 out <- compare.ZVrel(Vrel.gp$Jord, Vrel.gp$Teyah)
 
 summary(out)
 
## End(Not run)

Description

Subset (split) landmark coordinates via a grouping factor

Usage

coords.subset(A, group)

Arguments

Details

This function splits a set of landmark coordinates into subsets, as described by a factor. The result is a list of separate sets of landmarks.

Author(s)

Michael Collyer

Examples

## Not run: 
data(pupfish) 
group <- factor(paste(pupfish$Pop, pupfish$Sex))
levels(group)
new.coords <- coords.subset(A = pupfish$coords, group = group)
names(new.coords) # see the list levels
# group shape means
lapply(new.coords, mshape)

## End(Not run)

Description

An interactive function to define which landmarks should be linked to aid visualization.

Usage

define.links(spec, ptsize = 1, links = NULL)

Arguments

Details

Function takes a matrix of digitized landmark coordinates (e.g. from mshape) and allows the user to define pairs of landmarks to be linked, for visualization purposes. The output is a matrix to be used by plotAllSpecimens & plotRefToTarget option 'links='.

Note to RStudio users: if the 'zoom' of your RStudio window is set to something other than 100 sometimes the function can fail to identify the selected points ('warning: no point within 0.25 inches'). In these cases, set your zoom to 100

Selection

In the plot window select two landmarks that will be linked. In the console, the user will be prompted to continue linking landmarks to build the wireframe or end the session (typing y or n respectively). For 2D data in plot window, use LEFT mouse button to select landmarks. For 3D data in rgl window, use RIGHT mouse button (or command+LEFT for mac) to select landmarks.

Value

Function returns a matrix of which landmarks will be links to be used by plotAllSpecimens & plotRefToTarget option 'links='.

Author(s)

Emma Sherratt

See Also

plotAllSpecimens

plotRefToTarget

rgl-package (used in 3D plotting)

Description

An interactive function to define which landmarks should be assigned to each module (landmark partition).

Usage

define.modules(spec, nmodules)

Arguments

Details

Function takes a matrix of digitized landmark coordinates (e.g. from mshape) and allows the user to assign landmarks to each module. The output is a list of which landmarks belong in which partition, to be used by modularity.test or integration.test.

Selection in 2D

Choosing which landmarks will be included in each module involves landmark selection using a mouse in the plot window. The user is prompted to select each landmark in turn to be assigned to module 1: using the LEFT mouse button (or regular button for Mac users), click on the hollow circle to choose the landmark. Selected landmarks will be filled in. When all landmarks for module 1 are chosen, press 'esc', and then start selecting landmarks for module 2. Repeat until all modules are defined.

Selection in 3D

Choosing which landmarks will be included in each module involves landmark selection using a mouse in the rgl plot window. The user is prompted to select one or more landmarks. To do so, use the RIGHT mouse button (or command + LEFT button for Mac users), draw a rectangle around landmarks to select. Selected landmarks will be colored yellow. Then type into the console a letter (e.g. 1, 2, 3...) to assign selected landmark(s) to this module. Repeat until all landmarks are assigned to modules.

Value

Function returns a vector of which landmarks belong in which module (e.g. 1, 1, 1, 2, 2, 3, 3, 3, 2) to be used with modularity.test or integration.test.

Author(s)

Emma Sherratt

See Also

modularity.test and integration.test

rgl-package (used in 3D plotting)

Description

An interactive function to define which landmarks will "slide" along two-dimensional (2D) or three-dimensional (3D) curves.

Usage

define.sliders(landmarks, nsliders, surfsliders = NULL, write.file = TRUE)

Arguments

Details

Function takes a matrix of digitized landmark coordinates, such as made by digitize2d or digit.fixed, and helps user choose which landmarks will be treated as "sliders" in Generalized Procrustes analysis gpagen. This type of semilandmark "slides" along curves lacking known landmarks (see Bookstein 1997 for algorithm details). Each sliding semilandmark ("sliders") will slide between two designated points, along a line tangent to the specified curvature.

Defining landmarks is an interactive procedure (see below for 2D and 3D routines). The procedure is overlapping. For example: there are 5 landmarks (1:5), 1 and 5 are landmarks and 2,3,4 are sliders, the user must select '1' '2' '3', and then '2' '3' '4', and then '3' '4' '5'.

Selection in 2D

Choosing which landmarks will be sliders involves landmark selection using a mouse in the plot window. To define the sliders, for each sliding landmark along the curve in the format 'before-slider-after', using the LEFT mouse button (or regular button for Mac users), click on the hollow circle to choose the landmark in the following order:

1. Click to choose the first landmark between which semi-landmark will "slide",

2. Click to choose sliding landmark,

3. Click to choose the last landmark between which semi-landmark will "slide", Selected landmarks will be filled in and lines are drawn connecting the three landmarks, and will highlight the sliding semilandmark in red and the flanking landmarks in blue.

Selection in 3D

Choosing which landmarks will be sliders involves landmark selection using a mouse in the rgl plot window. With a standard 3-button (PC) buildtemplate uses:

1. the RIGHT mouse button (primary) to choose points to be defined as sliders,

2. the LEFT mouse button (secondary) is used to rotate mesh,

3. the mouse SCROLLER (third/middle) is used to zoom in and out.

NOTE: Digitizing functions on MACINTOSH computers using a standard 3-button mice works as specified. Macs using platform specific single button mice, XQuartz must be configured: go to Preferences > Input > tick "Emulate three button mouse":

1. press button to rotate 3D mesh,

2. press button while pressing COMMAND key to select points to be defined as sliders,

3. press button while pressing OPTION key to adjust mesh perspective.

4. the mouse SCROLLER or trackpad two finger scroll is used to zoom in an out.

To define the sliders, for each sliding landmark along the curve in the format 'before-slider-after':

1. Click to choose the first landmark between which semi-landmark will "slide",

2. Click to choose sliding landmark,

3. Click to choose the last landmark between which semi-landmark will "slide", Screen will show lines connecting the three landmarks, and will highlight the sliding semilandmark in red.

AUTO mode

The input 'landmarks' can be simply a vector of numbers corresponding to the "sliders" (semilandmarks) in the order they appear along a curve on the specimen. This can be made by c() or seq() or any other reasonable method.

If the sliders form a closed curve, then the function assumes that the first and last landmarks in the 'landmarks' vector are THE SAME are fixed (not sliders). e.g. if landmark 1 is a fixed landmark, and 2, 3 and 4 are semilandmarks, then sliders = c(1,2,3,4,1).

Value

Function returns a 'nsliders-x-3' matrix containing the landmark address of the curve sliders, indicating the landmarks between which the slider landmarks will "slide". If write.file = T the matrix is also written to working directory as "curveslide.csv". Matrix (or "curveslide.csv") is designed for use by gpagen during GPA.

Author(s)

Emma Sherratt & Dean Adams

References

Bookstein, F. J. 1997 Landmark Methods for Forms without Landmarks: Morphometrics of Group Differences in Outline Shape. Medical Image Analysis 1(3):225-243.

See Also

digitize2d, digit.fixed, gpagen, digit.curves

rgl-package (used in 3D plotting)

Examples

## (not run) Use interactive function in rgl window
 # data(scallops)
 # define.sliders(scallops$coorddata[,,1], nsliders=11,
 #   surfsliders = scallops$surfslide) 
 # here the first specimen is used for plotting purposes only
 
## Examples of AUTO mode 
 ## 1 curve of sliding semilandmark
 # Define sliders for scallopdata
 #sliders = define.sliders(c(5:16,1))

 ## 2 curves of sliding semilandmarks
 # Define sliders for 10 landmarks, where LMs 1, 5, and 10 fixed
 # 2, 3, and 4 are along a curve between 1 and 5
 # and 6, 7, 8, and 9 are along a curve between 5 and 10.
 #sliders = rbind(define.sliders(1:5), define.sliders(5:10))

Description

A function that "digitizes curves" by calculating equidistant two-dimensional or three-dimensional semilandmarks along a curve. These landmarks will be treated as "sliders" in Generalized Procrustes analysis gpagen. This type of semilandmark "slides" along curves lacking known landmarks (see Bookstein 1997 for algorithm details). Each sliding semilandmark ("sliders") will slide between two designated points, along a line tangent to the specified curvature, as specified by define.sliders.

Usage

digit.curves(start, curve, nPoints, closed = TRUE)

Arguments

Details

The function is based upon tpsDig2 'resample curve by length' for 2D data by James Rohlf (Rohlf 2015). The start of the curve is a fixed landmark on the curve that is equivalent (homologous) in each specimen in the sample (and will be treated as a fixed point during Procrustes Superimposition using gpagen). Then nPoints are calculated along the curve at equidistant points from the start to the end.

'curve' is a p-x-k matrix of 2D or 3D coordinates for a set of ordered points defining a curve. This can be the pixels of an outline calculated in ImageJ (save xy coordinates) or any other reasonable way of obtaining ordered coordinates along a curve (including sampling by hand using digit.fixed or digitize2d - but note that there should be more points defining the curve than nPoints in order to accurately calculate the semilandmarks).

If 'closed = T', the function returns the coordinates of the 'start' landmark plus nPoints. If 'closed = F', the function returns the coordinates of the 'start' landmark, plus nPoints and the end of the curve.

If unsure if the points defining the curve are ordered, then plot and color them using the rainbow function, e.g. plot(curve, pch=19, cex=0.1, col=rainbow(nrow(outline))), and it should be easy to visualize.

Value

Function returns a matrix of coordinates for nPoints equally spaced semilandmarks sampled along the curve (plus start and end if 'closed = F', or only including start if 'closed = T')

Author(s)

Emma Sherratt and Michael Collyer

References

Bookstein, F. J. 1997 Landmark Methods for Forms without Landmarks: Morphometrics of Group Differences in Outline Shape. Medical Image Analysis 1(3):225-243.

Rohlf, F.J., 2015. The tps series of software. Hystrix 26(1):9-12.

See Also

digit.fixed digitize2d

Description

An interactive function to digitize three-dimensional (3D) landmarks. Input for the function is either a matrix of vertex coordinates defining a 3D surface object or a mesh3d object as obtained from read.ply.

Usage

digit.fixed(spec, fixed, index = FALSE, ptsize = 1, center = TRUE)

Arguments

Details

Function for digitizing fixed three-dimensional landmarks. The user can later designate some as curve sliding semilandmarks, using the function define.sliders or through a semilandmark definition matrix.

To digitize 3D surface sliding semilandmarks the function digitsurface should be used instead.

For details on the full procedure for digitizing fixed 3D landmarks and surface sliding semilandmarks, see the relevant vignette by running vignette("geomorph.digitize3D").

NOTE: The function centers the mesh before digitizing by default (center = TRUE). If one chooses not to center, specimen may be difficult to manipulate in rgl window.

Digitizing

Digitizing is interactive. Once a point is selected, the user is asked if the system should keep or discard the selection (y/n). If "y", the user is asked to continue to select the next landmark. If "n" the removes the last chosen landmark, and the user is asked to select it again. This can be repeated until the user is comfortable with the landmark chosen.

To digitize with a standard 3-button mouse (PC):

1. the RIGHT mouse button (primary) to select points to be digitized,

2. the LEFT mouse button (secondary) is used to rotate mesh,

3. the mouse SCROLLER (third/middle) is used to zoom in and out.

NOTE: Digitizing functions on MACINTOSH computers using a standard 3-button mice works as specified. Macs using platform specific single button mice, XQuartz must be configured: go to Preferences > Input > tick "Emulate three button mouse":

1. press button to rotate 3D mesh,

2. press button while pressing COMMAND key to select vertex to be used as a landmark,

3. press button while pressing OPTION key to adjust mesh perspective.

4. the mouse SCROLLER or trackpad two finger scroll is used to zoom in an out.

NOTE: there is no pan (translate) functionality in rgl library for all platforms at this time.

Value

Function returns (if assigned to an object) and writes to the working directory an NTS file, containing the landmark coordinates. The file name corresponds to the name of the specimen. If index=FALSE function returns to the console an n x 3 matrix containing the x,y,z coordinates of the digitized landmarks. If index=TRUE, function returns a list:

Author(s)

Erik Otarola-Castillo & Emma Sherratt

See Also

read.ply

rgl-package (used in 3D plotting)

Description

An interactive function to digitize two-dimensional(2D) landmarks from .jpg files.

Usage

digitize2d(
  filelist,
  nlandmarks,
  scale = NULL,
  tpsfile,
  MultScale = FALSE,
  verbose = TRUE
)

Arguments

Details

This function may be used for digitizing 2D landmarks from jpeg images (.jpg). The user provides a list of image names, the number of landmarks to be digitized, and the name of an output TPS file. An option is included to allow the user to digitize a scale on each image to convert the landmark coordinates from pixels into meaningful units. Landmarks to be digitized can include both fixed landmarks and semi-landmarks, the latter of which are to be designated as "sliders" for subsequent analysis (see the function define.sliders).

The Digitizing Session

Digitizing landmarks from 2D photos requires that a scale bar is placed in the image in order to scale the coordinate data. The 'scale' option requires: a single number (e.g. 10) which means that the scale to be measured in all images is a 10mm scale bar; OR a vector the same length as the filelist containing a number for the scale of each image. If scale=NULL, then the digitized coordinates will not be scaled. This option is NOT recommended.

Users may digitize all specimens in one session, or may return at a later time to complete digitizing. In the latter case, the user provides the same filelist and TPS file and the function will determine where the user left off.

If specimens have missing landmarks, these can be incorporated during the digitizing process using the 'a' option as described below (a=absent).

Specimen Digitizing

Digitizing landmarks involves landmark selection using a mouse in the plot window, using the LEFT mouse button (or regular button for Mac users):

1. Digitize the scale bar (if requested) by selecting the two end points. Use a single click for start and end points. The user is asked whether the system should keep or discard the digitized scale bar.

2. Digitize each landmark with single click and the landmark is shown in red.

If verbose = TRUE, digitizing is interactive between landmark selection using a mouse and the R console. Once a landmark is selected, the user is asked if the system should keep or discard the selection (y/n/a). If "y", the user is asked to continue to select the next landmark. If "n", the user is asked to select it again.

To digitize a missing landmark, simply click on any location in the image. Then, when prompted to keep selection, choose 'a' (for absent). Missing landmarks can only be included during the digitizing process when verbose=TRUE.

If verbose = FALSE the digitizing of landmarks is continuous and uninterrupted. Here the user will not be prompted to approve each landmark selection.

At the end of digitizing, the landmark coordinates are written to a TPS file. By default, the x,y values are unscaled if a vector of scales is included, and the scale is returned on line SCALE= after each specimen x,y data. Optionally, one may have the coordinates pre-multiplied by scale by using the option MultScale=TRUE.

Value

Function returns a tps file containing the digitized landmark coordinates.

Author(s)

Dean Adams, Erik Otarola-Castillo and Emma Sherratt

See Also

readJPEG (for JPEG input)

Description

An interactive function to digitize three-dimensional (3D) landmarks on a surface lacking known landmarks. Input for the function is either a matrix of vertex coordinates defining a 3D surface object or a mesh3d object as obtained from read.ply.

Usage

digitsurface(spec, fixed, ptsize = 1, center = TRUE)

Arguments

Details

Function for digitizing fixed 3D landmarks and placing surface sliding semilandmarks using a previously created template. Following the selection of fixed points (see digitizing below), the function finds surface semilandmarks following the algorithm outlined in Gunz et al. (2005) and Mitteroecker and Gunz (2009). digitsurface finds the same number of surface semilandmarks as the template (created by buildtemplate) by downsampling the scanned mesh, after registering the template with the current specimen via GPA. A nearest neighbor algorithm is used to match template surface semilandmarks to mesh points of the current specimen. To use function digitsurface, the template must be constructed first, and 'template.txt' be in the working directory. Because template matching is based on the correspondence of fixed landmark points in the template and the specimen, a minimum of four fixed landmarks must be used.

For details on the full procedure for digitizing fixed 3D landmarks and surface sliding semilandmarks, see the relevant vignette by running vignette("geomorph.digitize3D").

NOTE: Function centers the mesh before digitizing by default (center = TRUE). If one chooses not to center, specimen may be difficult to manipulate in rgl window.

Digitizing

Digitizing of fixed landmarks is interactive. Once a point is selected, the user is asked if the system should keep or discard the selection (y/n). If "y", the user is asked to continue to select the next landmark. If "n" the removes the last chosen landmark, and the user is asked to select it again. This can be repeated until the user is comfortable with the landmark chosen.

To digitize with a standard 3-button (PC):

1. the RIGHT mouse button (primary) to select points to be digitized,

2. the LEFT mouse button (secondary) is used to rotate mesh,

3. the mouse SCROLLER (third/middle) is used to zoom in and out.

NOTE: Digitizing functions on MACINTOSH computers using a standard 3-button mice works as specified. Macs using platform specific single button mice, XQuartz must be configured: go to Preferences > Input > tick "Emulate three button mouse":

1. press button to rotate 3D mesh,

2. press button while pressing COMMAND key to select vertex to be used as a landmark,

3. press button while pressing OPTION key to adjust mesh perspective.

4. the mouse SCROLLER or trackpad two finger scroll is used to zoom in an out.

NOTE: there is no pan (translate) functionality in rgl library for all platforms at this time.

AUTO mode

The function as described above (for interactive mode) calls digit.fixed, prompting the user to select fixed landmarks in the rgl window. However if the user has digitized these fixed landmark elsewhere (e.g., in other software), then the input for parameter 'fixed' can be a p-x-k matrix of 3D coordinates. In this case, the function the function will automatically use these landmarks and fit the template of sliding semilandmarks.

Value

Function returns (if assigned to an object) and writes to the working directory an NTS file, containing the landmark coordinates. The file name corresponds to the name of the specimen.

Author(s)

Erik Otarola-Castillo & Emma Sherratt

References

Gunz P, Mitteroecker P, & Bookstein FJ (2005) Semilandmarks in Three Dimensions. Modern Morphometrics in Physical Anthropology, ed Slice DE (Springer-Verlag, New York), pp 73-98.

Mitteroecker P & Gunz P (2009) Advances in Geometric Morphometrics. Evolutionary Biology 36(2):235-247.

See Also

buildtemplate

read.ply

digit.fixed

rgl-package (used in 3D plotting)

Description

An interactive function to remove landmarks from a 3D template file.

Usage

editTemplate(template, fixed, n)

Arguments

Details

The function edits a 'template.txt' file made by buildtemplate, which must be in the current working directory, and which is overwritten. Use read.table("template.txt", header = T) to read in the template first.

Selection

Choosing which landmarks will be deleted involves landmark selection using a mouse in the rgl plot window. With a standard 3-button (PC) buildtemplate uses:

1. the RIGHT mouse button (primary) to choose points to be deleted,

2. the LEFT mouse button (secondary) is used to rotate mesh,

3. the mouse SCROLLER (third/middle) is used to zoom in and out.

NOTE: Digitizing functions on MACINTOSH computers using a standard 3-button mice works as specified. Macs using platform specific single button mice, XQuartz must be configured: go to Preferences > Input > tick "Emulate three button mouse":

1. press button to rotate 3D mesh,

2. press button while pressing COMMAND key to select points to be deleted,

3. press button while pressing OPTION key to adjust mesh perspective.

4. the mouse SCROLLER or trackpad two finger scroll is used to zoom in an out.

Value

Function returns a matrix containing the x,y,z coordinates of the new template landmarks. Function also writes to the working directory 'template.txt' containing the x,y,z coordinates of the updated template

Author(s)

Erik Otarola-Castillo & Emma Sherratt

See Also

rgl-package (used in 3D plotting)

Description

A function for estimating the locations of missing landmarks

Usage

estimate.missing(A, method = c("TPS", "Reg"))

Arguments

Details

The function estimates the locations of missing landmarks for incomplete specimens in a set of landmark configurations, where missing landmarks in the incomplete specimens are designated by NA in place of the x,y,z coordinates. Two distinct approaches are implemented.

The first approach (method = "TPS") uses the thin-plate spline to interpolate landmarks on a reference specimen to estimate the locations of missing landmarks on a target specimen. Here, a reference specimen is obtained from the set of specimens for which all landmarks are present, Next, each incomplete specimen is aligned to the reference using the set of landmarks common to both. Finally, the thin-plate spline is used to estimate the locations of the missing landmarks in the target specimen (Gunz et al. 2009).

The second approach (method = "Reg") is multivariate regression. Here each landmark with missing values is regressed on all other landmarks for the set of complete specimens, and the missing landmark values are then predicted by this linear regression model. Because the number of variables can exceed the number of specimens, the regression is implemented on scores along the first set of PLS axes for the complete and incomplete blocks of landmarks (see Gunz et al. 2009). Note, however, that a minimum of k*m+k specimens are required to estimate m missing landmarks (of k-dimension) in any one specimen using the regression method. More generally, if the number of missing landmarks approaches the number of reference specimens used to estimate them, estimation will become increasingly imprecise with the regression method. Additionally, the location of missing landmarks (contiguous versus disparate in location) can also influence the precision of estimation. The user should be aware that the function will produce results but the results from the regression method might be influenced by the number of specimens, the number of total landmarks, and the number and location of missing landmarks in any one specimen. It might be wise to compare multiple methods for specific cases, if uncertain about the precision of estimation.

One can also exploit bilateral symmetry to estimate the locations of missing landmarks. Several possibilities exist for implementing this approach (see Gunz et al. 2009). Example R code for one implementation is found in Claude (2008).

NOTE: Because all geometric morphometric analyses and plotting functions implemented in geomorph require a full complement of landmark coordinates, the alternative to estimating the missing landmark coordinates is to proceed with subsequent analyses EXCLUDING specimens with missing values.

Value

Function returns an array (p x k x n) of the same dimensions as input A, including coordinates for the target specimens (the original landmarks plus the estimated coordinates for the missing landmarks). If the input is a geomorphShapes object, this is returned with the original and estimated coordinates in $landmarks In both cases, these data need to be Procrustes Superimposed prior to analysis, including sliding of semilandmarks (see gpagen).

Author(s)

Dean Adams

References

Claude, J. 2008. Morphometrics with R. Springer, New York.

Bookstein, F. L., K. Schafer, H. Prossinger, H. Seidler, M. Fieder, G. Stringer, G. W. Weber, J.-L. Arsuaga, D. E. Slice, F. J. Rohlf, W. Recheis, A. J. Mariam, and L. F. Marcus. 1999. Comparing frontal cranial profiles in archaic and modern Homo by morphometric analysis. Anat. Rec. (New Anat.) 257:217-224.

Gunz, P., P. Mitteroecker, S. Neubauer, G. W. Weber, and F. L. Bookstein. 2009. Principles for the virtual reconstruction of hominin crania. J. Hum. Evol. 57:48-62.

Examples

## Not run: 
data(plethodon)
plethland <- plethodon$land
  plethland[3,,2] <- plethland[8,,2] <- NA  #create missing landmarks
  plethland[3,,5] <- plethland[8,,5] <- plethland[9,,5] <- NA  
  plethland[3,,10] <- NA  
  
estimate.missing(plethland, method = "TPS")
estimate.missing(plethland, method = "Reg")

## End(Not run)

Description

Function performs extended PGLS that enables more than one individual per species to be included in a phylogenetic least squares analysis.

Usage

extended.pgls(
  f1,
  phy = NULL,
  Cov = NULL,
  species = NULL,
  delta = 0.001,
  gamma = c("sample", "equal"),
  iter = 999,
  seed = NULL,
  int.first = FALSE,
  turbo = FALSE,
  Parallel = FALSE,
  verbose = FALSE,
  data = NULL,
  print.progress = TRUE,
  ...
)

Arguments

Details

The function performs linear models in a phylogenetic context in a manner that can accommodate multiple individuals per species and can accommodate high-dimensional datasets. The approach utilizes a hierarchical linear model, an expanded phylogenetic covariance matrix, and permutation procedures to obtain empirical sampling distributions and effect sizes for model effects that can evaluate differences in intraspecific trends across species for both univariate and multivariate data, while conditioning them on the phylogeny (Adams and Collyer 2024).

Data input is specified by several components, including minimally: 1) a formula (e.g., y ~ X), where 'y' specifies the response variables (shape data), and 'X' contains one or more independent variables (discrete or continuous). The response matrix 'Y' can be either in the form of a two-dimensional data matrix of dimension (n x [p x k]), or a 3D array (p x n x k). It is assumed that the landmarks have previously been aligned using Generalized Procrustes Analysis (GPA) [e.g., with gpagen].

2) An argument for species is also required, which specifies the species to which each observation belongs.

3) Either phy or COV must be included, which specifies the expected phylogenetic non-independence among species. Here, 'phy' represents a phylogenetic tree of class 'phylo'. In this case, the expected covariances among species will be estimated under a Brownian motion model of evolution. Alternatively, 'COV' is a phylogenetic covariance matrix obtained previously by the user, and may represent the expected covariances among species under any model of evolutionary change (BM, OU, etc.).

From the phylogenenetic covariance matrix, a phylogenetic transformation matrix is obtained and used to transform the Y variables and the model design matrices constructed from X, and parameter estimates are obtained. Next, variance components for all model terms are calculated, using a combination of type II and type III sums of squares, so that within-species effects may be isolated from interspecies effects (see Adams and Collyer 2024). Likewise, residuals from the model are permuted in such a manner as to isolate intra- and inter-species effects, to account for fact that multiple observations are included for each species. Additional statistical and philosophical details are found in Adams and Collyer 2024.

Value

An object of class lm.rrpp.ws is a list containing the following

Author(s)

Dean Adams

References

Adams, D.C and M.L Collyer. 2024. Extending phylogenetic regression models for comparing within-species patterns across the tree of life. Methods in Ecology and Evolution. DOI: 10.1111/2041-210X.14438

See Also

lm.rrpp.ws

Examples

## Not run: 
data(pupfish.ws)

# With phylogeny
fit <- extended.pgls(f1 = coords~Species * Sex + Population, 
  data = pupfish.ws, species = "Species",
  phy = pupfish.ws$phy) 
  
anova(fit) 
 
# multivariate stats 
fit.mult <- manova.update(fit, PC.no = 40)
summary(fit.mult, test = "Wilks") 

# With phylogenetic covariance matrix

fit2 <- extended.pgls(f1 = coords ~ Species * Sex + Population, 
  data = pupfish.ws, species = "Species",
  Cov = pupfish.ws$Cov)
  
anova(fit2)

#sultivariate stats
fit2.mult <- manova.update(fit2, PC.no = 40)
summary(fit2.mult, test = "Wilks") 

## End(Not run)

Description

A function to identify which specimen lies closest to the estimated mean shape for a set of Procrustes shape variables.

Usage

findMeanSpec(A)

Arguments

Details

Function takes a 3D array of Procrustes shape variables (such as made by gpagen, calculates the distance of each to the estimated mean shape, and returns the name and address of the closest specimen. This function can be used to identify the specimen to be used by warpRefMesh.

Value

Function returns the name and address of the specimen closest to the mean of the set of Procrustes shape variables.

Author(s)

Emma Sherratt

See Also

warpRefMesh

Description

A function for rotating a subset of landmarks so that the articulation angle between subsets is constant

Usage

fixed.angle(
  A,
  art.pt = NULL,
  angle.pts.1,
  angle.pts.2,
  rot.pts = NULL,
  angle = 0,
  degrees = FALSE
)

Arguments

Details

This function standardizes the angle between two subsets of landmarks for a set of specimens. The approach assumes a simple hinge-point articulation between the two subsets, and rotates all specimens such that the angle between landmark subsets is equal across specimens (see Adams 1999). As a default, the mean angle is used, though the user may specify an additional amount by which this may be augmented. To quantify the angle, users may specify a single landmark in each subset as angle endpoints, or may specify a set of landmarks. If the latter, the centroid of those points is used.

Presently, the function is only implemented for two-dimensional landmark data.

Value

Function returns a (p x k x n) array of landmark coordinates.

Author(s)

Dean Adams and Michael Collyer

References

Adams, D. C. 1999. Methods for shape analysis of landmark data from articulated structures. Evolutionary Ecology Research. 1:959-970.

Examples

## Not run: 
#Example using Plethodon
#Articulation point is landmark 1, rotate mandibular landmarks (2-5) 
# relative to cranium

data(plethspecies) 
# Using specific points:
newLM1 <- fixed.angle(plethspecies$land,
art.pt = 1, angle.pts.1 = 5, 
angle.pts.2 = 6, rot.pts = c(2,3,4,5))
Y.gpa1 <- gpagen(newLM1)
plot(Y.gpa1, mean = FALSE)

# Using centroids from subsets
newLM2 <- fixed.angle(plethspecies$land, art.pt = 1, 
angle.pts.1 = c(1, 6:11), 
angle.pts.2 = 2:5, 
rot.pts = NULL, angle = 20, 
degrees = TRUE) # rotated points same as second partition
Y.gpa2 <- gpagen(newLM2)
plot(Y.gpa2, mean = FALSE)

## End(Not run)

Description

A list similar to a data frame to facilitate analysis of shape data.

Usage

geomorph.data.frame(...)

Arguments

Details

This function produces a list that can be used like a data frame in other analytical functions. The purpose is similar to the function, data.frame, but without the constraint that data must conform to an n (observations) x p (variables) matrix. Rather, the list produced is constrained only by n. List objects can be Procrustes shape variables, matrices, variables, distance matrices, and phylogenetic trees. Results from gpagen can be directly imported into a geomorph.data.frame to utilize the coordinates and centroid size as variables. (See Examples)

Author(s)

Michael Collyer

Examples

## Not run: 
data(plethodon) 
Y.gpa <- gpagen(plethodon$land, PrinAxes = FALSE)
gdf <- geomorph.data.frame(Y.gpa)
attributes(gdf)

gdf <- geomorph.data.frame(Y.gpa, species = plethodon$species, 
site = plethodon$site)
attributes(gdf)

# Using geomorph.data.frame to facilitate analysis
anova(procD.lm(coords ~ Csize + species * site, data = gdf))

## End(Not run)

Description

Function quantifies the overall level of morphological integration for a set of Procrustes shape variables

Usage

globalIntegration(A, ShowPlot = TRUE)

Arguments

Details

The function quantifies the overall level of morphological integration for a set of Procrustes shape coordinates. It is assumed that the landmarks have previously been aligned using Generalized Procrustes Analysis (GPA) [e.g., with gpagen]. Based on the set of aligned specimens, the function estimates the set of bending energies at various spatial scales, and plots the log of the variance of the partial warps versus the log of their corresponding bending energies (Bookstein 2015). The slope of a regression of these data provides information regarding the degree of overall morphological integration (or lack thereof).

A slope of negative one corresponds to self-similarity, implying that patterns of shape variation are similar across spatial scales. Steeper slopes (i.e., those more extreme than -1.0) correspond to data that are globally integrated, while shallower slopes (between -1 and 0) correspond to data that are 'disintegrated (see Bookstein 2015). Isotropic data will have an expected slope of zero.

Author(s)

Dean Adams

References

Bookstein, F. L. 2015. Integration, disintegration, and self-similarity: Characterizing the scales of shape variation in landmark data. Evol. Biol.42(4): 395-426.

Examples

## Not run: 
data(plethodon) 
Y.gpa <- gpagen(plethodon$land)    #GPA-alignment    

globalIntegration(Y.gpa$coords)

## End(Not run)

Description

Evaluation of measurement error for two or more multivariate measurements, for common research subjects.

Usage

gm.measurement.error(
  coords,
  subjects,
  replicates,
  groups = NULL,
  data,
  iter = 999,
  seed = NULL,
  multivariate = FALSE,
  use.PCs = TRUE,
  tol = 0.001,
  Parallel = FALSE,
  turbo = TRUE,
  print.progress = FALSE,
  verbose = FALSE,
  ...
)

Arguments

Details

Function performs analyses concerned with the repeatability (reliability) of multivariate data (measurements) collected from the same research subjects. Although there is no requirement for repeated measurements on all research subjects, the analysis assumes that multiple observations are made.

This function performs analyses as described in Collyer and Adams (in press) to assess systematic and random components of measurement error (ME). It basically performs ANOVA with RRPP, but with different restricted randomization strategies. The reliability of research subject variation can be considered by restricting randomization within replicates; the consistency of replicate measures can be considered by restricting randomization within subjects. Inter-subject variation remains constant across all random permutations within subjects and inter-replicate variation remains constant across all random permutations within replicates. Type II sums of squares and cross-products (SSCP) are calculated to assure conditional estimation.

The results include univariate-like statistics based on dispersion of values and eigenanalysis performed on a signal to noise matrix product of SSCP matrices (sensu Bookstein and Mitteroecker, 2014) including the inverse of the random component of ME and the systematic component of ME. The multivariate test is a form of multivariate ANOVA (MANOVA), using RRPP to generate sampling distributions of the major eigenvalue (Roy's maximum root).

Value

Objects of class "measurement.error" return the same objects as a lm.rrpp fit, plus a list of the following:

Author(s)

Michael Collyer and Dean Adams

References

Collyer, M.L. and D.C. Adams. 2024. Interrogating Random and Systematic Measurement Error in Morphometric Data. Evolutionary Biology.

Bookstein, F.L., & Mitteroecker, P. (2014). Comparing covariance matrices by relative eigenanalysis, with applications to organismal biology. Evolutionary biology, 41(2), 336-350.

See Also

lm.rrpp.ws, manova.update

Examples

## Not run: 
# Measurement error analysis on simulated data of fish shapes

data(fishy)
fishy$coordsarray <- arrayspecs(fishy$coords, p = 11, k = 2)  #make 3D array

# Example two digitization replicates of the same research subjects
rep1 <- matrix(fishy$coords[1,], 11, 2, byrow = TRUE)
rep2 <- matrix(fishy$coords[61,], 11, 2, byrow = TRUE)
plot(rep1, pch = 16, col = gray(0.5, alpha = 0.5), cex = 2, asp = 1)
points(rep2, pch = 16, col = gray(0.2, alpha = 0.5), cex = 2, asp = 1)

# Analysis unconcerned with groups 

ME1 <- gm.measurement.error(
  coords = "coordsarray",
  subjects = "subj",
  replicates = "reps",
  data = fishy)

anova(ME1)
ICCstats(ME1, subjects = "Subjects", with_in = "Systematic ME")
plot(ME1)

# Analysis concerned with groups 

ME2 <- gm.measurement.error(
  coords = "coordsarray",
  subjects = "subj",
  replicates = "reps",
  groups = "groups",
  data = fishy)
  
anova(ME2)
ICCstats(ME2, subjects = "Subjects", 
  with_in = "Systematic ME", groups = "groups")
P <- plot(ME2)
focusMEonSubjects(P, subjects = 18:20, shadow = TRUE)

#heat map of inter-subject variability
int.var <- interSubVar(ME2, type = "var")
plot(int.var)

## End(Not run)

Description

Function performs principal components analysis (PCA) or phylogenetically-aligned components (PaCA) on Procrustes shape coordinates.

Usage

gm.prcomp(
  A,
  phy = NULL,
  align.to.phy = FALSE,
  unit.tree = TRUE,
  GLS = FALSE,
  transform = FALSE,
  ...
)

Arguments

Details

The function performs a series of ordinations, taking into account, phylogeny, if desired. There are two main types of ordinations: principal components analysis (PCA) and phylogenetically- aligned components analysis (PaCA). Both of these have two variants: centering and projection via ordinary least-squares (OLS) or via generalized least-squares (GLS). The name, "gm.prcomp", references that this function performs much like prcomp, in terms of arguments and output, but this function is quite a bit more diverse. This function has the capability of performing analyses generally referred to as:

PCA

Standard PCA based on OLS-centering and projection of data.

Phylomorphospace

Standard PCA with estimated ancestral states and phylogenetic branches projected into ordination plots.

phyloPCA

PCA based on GLS-centering and projection of data. Also possible to project ancestral states into plots. Note that if transformed GLS-residuals are used for projection, the ancestral states might not appear logical, as the projection is independent of phylogeny. With OLS-centering, a phyloPCA as described by Revell (2009) is produced.

PaCA

Phylogenetically-aligned component analysis. Data are aligned to an axis of greatest phylogenetic signal rather than axis of greatest dispersion. This analysis can use either OLS- or GLS-centering and projection. Phylogenetic signal is strongest in the first few components of the OLS approach. This analysis will make little sense with GLS-centering and projection of transformed residuals, since phylogenetic signal is removed the transformed data. See Collyer and Adams (2021) for more details. For greater flexibility for type of residuals and projection of trees, use ordinate. See Collyer and Adams (2021) for details.

phy

Whether a phylogeny and estimated ancestral states are considered in plots.

align.to.phy

Whether components are aligned to phylogenetic signal (rather than principal axes).

GLS

Whether to use GLS-centering and estimation of covariance matrix.

transform

Whether to transform GLS-residuals (making them independent of phylogeny and an orthogonal projection from the transformed data space, as opposed to an oblique projection from the untransformed data space).

PLOTTING: Contrary to previous geomorph implementations, gm.prcomp does not produce plots. For plotting gm.prcomp class objects combine plot.gm.prcomp and picknplot.shape following the examples below.

SUMMARY STATISTICS: For principal component plots, the traditional statistics to summarize the analysis include eigenvalues (variance by component), proportion of variance by component, and cumulative proportion of variance. When data are aligned to a phylogenetic covariance matrix, the statistics are less straightforward. A summary of of such an analysis (performed with summary.gm.prcomp) will produce these additional statistics:

Singular Value

Rather than eigenvalues, the singular values from singular value decomposition of the cross-product of the scaled phylogenetic covariance matrix and the data.

Proportion of Covariance

Each component's singular value divided by the sum of singular values. The cumulative proportion is also returned. Note that these values do not explain the amount of covariance between phylogeny and data, but explain the distribution of the covariance. Large proportions can be misleading.

RV by Component

The partial RV statistic by component. Cumulative values are also returned. The sum of partial RVs is Escoffier's RV statistic, which measures the amount of covariation between phylogeny and data. Caution should be used in interpreting these values, which can vary with the number of observations and number of variables. However, the RV is more reliable than proportion of singular value for interpretation of the strength of linear association for phylogenetically-aligned components.

Tips or Ancestors Dispersion

The variances of points by component for tip data and estimated ancestral character states, after projection. These values will differ from variances from PCA with GLS estimation, as the "Importance of Components" weights variances by phylogenetic covariances. Dispersion statistics correspond to the amount of scatter in plots of component scores.

NOTE: The plot.gm.prcomp function performs the same plotting that was previously possible with plotTangentSpace and plotGMPhyloMorphoSpace, which have now been deprecated.

Value

An object of class "gm.prcomp" contains a list of results for each of the PCA approaches implemented. Each of these lists includes the following components:

Author(s)

Antigoni Kaliontzopoulou, Michael Collyer, & Dean Adams

References

Collyer, M.L and D.C. Adams, 2021. Phylogenetically Aligned component analysis. Methods in Ecology and Evolution, 12: 369-372.

Revell, L. J. (2009). Size-correction and principal components for interspecific comparative studies. Evolution, 63: 3258-3268.

See Also

plot.gm.prcomp

picknplot.shape

ordinate

Examples

## Not run: 
 data(plethspecies) 
 Y.gpa <- gpagen(plethspecies$land)    #GPA-alignment
 
 ###  Traditional PCA 
 PCA <- gm.prcomp(Y.gpa$coords)
 summary(PCA)
 plot(PCA, main = "PCA")
 plot(PCA, main = "PCA", flip = 1) # flip the first axis
 plot(PCA, main = "PCA", axis1 = 3, axis2 = 4) # change PCs viewed
 
 ### Phylomorphospace - PCA with phylogeny (result is same as above, 
 ### but with estimated ancestral states projected into plot)
 PCA.w.phylo <- gm.prcomp(Y.gpa$coords, phy = plethspecies$phy)
 summary(PCA.w.phylo)
 plot(PCA.w.phylo, phylo = TRUE, main = "PCA.w.phylo")
 
 ### Phylogenetic PCA - PCA based on GLS-centering and projection
 # This is the same as the method described by Revell (2009)
 phylo.PCA <- gm.prcomp(Y.gpa$coords, phy = plethspecies$phy, GLS = TRUE)
 summary(phylo.PCA)
 plot(phylo.PCA, phylo = TRUE, main = "phylo PCA")
 
 ### Phylogenetic PCA - PCA based on GLS-centering and transformed 
 # projection
 # This produces a PCA independent of phylogeny
 phylo.tPCA <- gm.prcomp(Y.gpa$coords, phy = plethspecies$phy, 
 GLS = TRUE, transform = TRUE)
 summary(phylo.tPCA)
 plot(phylo.tPCA, phylo = TRUE, main = "phylo PCA")
 
 ### PaCA - Alignment of data to physlogenetic signal rather than axis of 
 ### greatest variation, like in PCA
 
 # OLS method (rotation of PCA)
 PaCA.ols <- gm.prcomp(Y.gpa$coords, phy = plethspecies$phy, 
   align.to.phy = TRUE)
 summary(PaCA.ols)
 plot(PaCA.ols, phylo = TRUE, main = "PaCA using OLS")
 
 # GLS method (rotation of Phylogenetic PCA)
 PaCA.gls <- gm.prcomp(Y.gpa$coords, phy = plethspecies$phy, 
 align.to.phy = TRUE, GLS = TRUE)
 summary(PaCA.gls)
 plot(PaCA.gls, phylo = TRUE, main = "PaCA using GLS")
 
 # GLS method (rotation of Phylogenetic PCA with transformed data)
 PaCA.gls <- gm.prcomp(Y.gpa$coords, phy = plethspecies$phy, 
 align.to.phy = TRUE, GLS = TRUE, transform = TRUE)
 summary(PaCA.gls)
 plot(PaCA.gls, phylo = TRUE, 
   main = "PaCA using GLS and transformed projection")
 
 
 ### Advanced Plotting
 gps <- as.factor(c(rep("gp1", 5), rep("gp2", 4))) # Two random groups
 par(mar=c(2, 2, 2, 2))
 plot(PaCA.ols, pch=22, cex = 1.5, bg = gps, phylo = TRUE) 
 # Modify options as desired
 #  Add things as desired using standard R plotting
 text(par()$usr[1], 0.1*par()$usr[3], labels = "PC1 - 45.64%", 
   pos = 4, font = 2)
 text(0, 0.95*par()$usr[4], labels = "PC2 - 18.80%", pos = 4, font = 2)
 legend("topleft", pch=22, pt.bg = unique(gps), legend = levels(gps))
 
 ## 3D plot with a phylogeny and time on the z-axis
  plot(PCA.w.phylo, time.plot = TRUE)
  plot(PCA.w.phylo, time.plot = TRUE, bg = "red", 
     phylo.par = list(tip.labels = TRUE, 
  tip.txt.cex = 2, edge.color = "blue", edge.width = 2))
 
## End(Not run)

Description

A general function to perform Procrustes analysis of two- or three-dimensional landmark data that can include both fixed landmarks and sliding semilandmarks

Usage

gpagen(
  A,
  curves = NULL,
  surfaces = NULL,
  rot.pts = NULL,
  PrinAxes = TRUE,
  max.iter = NULL,
  tol = 1e-04,
  ProcD = FALSE,
  approxBE = FALSE,
  sen = 0.5,
  Proj = TRUE,
  verbose = FALSE,
  print.progress = TRUE,
  Parallel = FALSE
)

Arguments

Details

The function performs a Generalized Procrustes Analysis (GPA) on two-dimensional or three-dimensional landmark coordinates. The analysis can be performed on fixed landmark points, semilandmarks on curves, semilandmarks on surfaces, or any combination. If data are provided in the form of a 3D array, all landmarks and semilandmarks are contained in this object. If this is the only component provided, the function will treat all points as if they were fixed landmarks. To designate some points as semilandmarks, one uses the "curves =" or "surfaces =" options (or both). To include semilandmarks on curves, a matrix defining which landmarks are to be treated as semilandmarks is provided using the "curves =" option. This matrix contains three columns that specify the semilandmarks and two neighboring landmarks which are used to specify the tangent direction for sliding. The matrix may be generated using the function define.sliders). Likewise, to include semilandmarks on surfaces, one must specify a vector listing which landmarks are to be treated as surface semilandmarks using the "surfaces=" option. The "ProcD = FALSE" option (the default) will slide the semilandmarks based on minimizing bending energy, while "ProcD = TRUE" will slide the semilandmarks along their tangent directions using the Procrustes distance criterion. The Procrustes aligned specimens may be projected into tangent space using the "Proj = TRUE" option. The function also outputs a matrix of pairwise Procrustes Distances, which correspond to Euclidean distances between specimens in tangent space if "Proj=TRUE", or to the geodesic distances in shape space if "Proj=FALSE". NOTE: Large datasets may exceed the memory limitations of R.

Generalized Procrustes Analysis (GPA: Gower 1975, Rohlf and Slice 1990) is the primary means by which shape variables are obtained from landmark data (for a general overview of geometric morphometrics see Bookstein 1991, Rohlf and Marcus 1993, Adams et al. 2004, Zelditch et al. 2012, Mitteroecker and Gunz 2009, Adams et al. 2013). GPA translates all specimens to the origin, scales them to unit-centroid size, and optimally rotates them (using a least-squares criterion) until the coordinates of corresponding points align as closely as possible. The resulting aligned Procrustes coordinates represent the shape of each specimen, and are found in a curved space related to Kendall's shape space (Kendall 1984). Typically, these are projected into a linear tangent space yielding Kendall's tangent space coordinates (i.e., Procrustes shape variables), which are used for subsequent multivariate analyses (Dryden and Mardia 1993, Rohlf 1999). Additionally, any semilandmarks on curves and surfaces are slid along their tangent directions or tangent planes during the superimposition (see Bookstein 1997; Gunz et al. 2005). Presently, two implementations are possible: 1) the locations of semilandmarks can be optimized by minimizing the bending energy between the reference and target specimen (Bookstein 1997), or by minimizing the Procrustes distance between the two (Rohlf 2010). Note that specimens are NOT automatically reflected to improve the GPA-alignment.

The generic functions, print, summary, and plot all work with gpagen. The generic function, plot, calls plotAllSpecimens.

Notes for geomorph 4.0

Starting with geomorph version 4.0, it is possible to use approximated thin-plate spline (TPS) mapping to estimate bending energy (when bending energy is used for sliding semi-landmarks). Approximated TPS has some computational benefit for large data sets (more GPA iterations are possible in a shorter time), but one should make sure first that results are reasonable. Approximated TPS uses a subset of points (semilandmarks) to estimate bending energy locally, where points are free to slide. Some subsets might be misrepresentative of the global (full configuration) bending energy, and strange results are possible. A comparison between TPS and approximated TPS outcomes could be tried with a few specimens to verify consistency of results before using approximated TPS on a large data set.

Choosing to use approximated TPS with only a few sliders is dangerous, as the subset of points might be too small to be reliable. Approximated TPS will work best for data sets with a sufficient number of surface landmarks. If a data set is nearly 100 % semilandmarks, TPS and approximated TPS should converge, but no time savings is likely.

Notes for geomorph 3.0

Compared to older versions of geomorph, users might notice subtle differences in Procrustes shape variables when using semilandmarks (curves or surfaces). This difference is a result of using recursive updates of the consensus configuration with the sliding algorithms (minimized bending energy or Procrustes distances). (Previous versions used a single consensus through the sliding algorithms.) Shape differences using the recursive updates of the consensus configuration should be highly correlated with shape differences using a single consensus during the sliding algorithm, but rotational "flutter" can be expected. This should have no qualitative effect on inferential analyses using Procrustes residuals.

Value

An object of class gpagen returns a list with the following components:

Author(s)

Dean Adams and Michael Collyer

References

Adams, D. C., F. J. Rohlf, and D. E. Slice. 2004. Geometric morphometrics: ten years of progress following the 'revolution'. It. J. Zool. 71:5-16.

Adams, D. C., F. J. Rohlf, and D. E. Slice. 2013. A field comes of age: Geometric morphometrics in the 21st century. Hystrix.24:7-14.

Bookstein, F. L. 1991. Morphometric tools for landmark data: Geometry and Biology. Cambridge Univ. Press, New York.

Bookstein, F. L. 1997. Landmark methods for forms without landmarks: morphometrics of group differences in outline shape. 1:225-243.

Dryden, I. L., and K. V. Mardia. 1993. Multivariate shape analysis. Sankhya 55:460-480.

Gower, J. C. 1975. Generalized Procrustes analysis. Psychometrika 40:33-51.

Gunz, P., P. Mitteroecker, and F. L. Bookstein. 2005. semilandmarks in three dimensions. Pp. 73-98 in D. E. Slice, ed. Modern morphometrics in physical anthropology. Klewer Academic/Plenum, New York.

Kendall, D. G. 1984. Shape-manifolds, Procrustean metrics and complex projective spaces. Bulletin of the London Mathematical Society 16:81-121.

Mitteroecker, P., and P. Gunz. 2009. Advances in geometric morphometrics. Evol. Biol. 36:235-247.

Rohlf, F. J., and D. E. Slice. 1990. Extensions of the Procrustes method for the optimal superimposition of landmarks. Syst. Zool. 39:40-59.

Rohlf, F. J., and L. F. Marcus. 1993. A revolution in morphometrics. Trends Ecol. Evol. 8:129-132.

Rohlf, F. J. 1999. Shape statistics: Procrustes superimpositions and tangent spaces. Journal of Classification 16:197-223.

Rohlf, F. J. 2010. tpsRelw: Relative warps analysis. Version 1.49. Department of Ecology and Evolution, State University of New York at Stony Brook, Stony Brook, NY.

Zelditch, M. L., D. L. Swiderski, H. D. Sheets, and W. L. Fink. 2012. Geometric morphometrics for biologists: a primer. 2nd edition. Elsevier/Academic Press, Amsterdam.

Examples

## Not run: 
# Example 1: fixed points only
data(plethodon) 
Y.gpa <- gpagen(plethodon$land, PrinAxes = FALSE)
summary(Y.gpa)
plot(Y.gpa)

# Example 2: points and semilandmarks on curves
data(hummingbirds)

###Slider matrix
hummingbirds$curvepts

# Using bending energy for sliding
Y.gpa <- gpagen(hummingbirds$land, curves = hummingbirds$curvepts,
ProcD = FALSE)   
summary(Y.gpa)
plot(Y.gpa)


# Using Procrustes Distance for sliding
Y.gpa <- gpagen(hummingbirds$land, curves = hummingbirds$curvepts,
ProcD = TRUE)   
summary(Y.gpa)
plot(Y.gpa)


# Example 3: points, curves and surfaces
data(scallops)

# Using Procrustes Distance for sliding
Y.gpa <- gpagen(A = scallops$coorddata, curves = scallops$curvslide, 
surfaces = scallops$surfslide)
# NOTE can summarize as: summary(Y.gpa)
# NOTE can plot as: plot(Y.gpa) 

## End(Not run)

Description

Function produces a list of parameter changes within plotRefToTarget.

Usage

gridPar(
  pt.bg = "gray",
  pt.size = 1.5,
  link.col = "gray",
  link.lwd = 2,
  link.lty = 1,
  out.col = "gray",
  out.cex = 0.1,
  tar.pt.bg = "black",
  tar.pt.size = 1,
  tar.link.col = "black",
  tar.link.lwd = 2,
  tar.link.lty = 1,
  tar.out.col = "black",
  tar.out.cex = 0.1,
  n.col.cell = 20,
  grid.col = "black",
  grid.lwd = 1,
  grid.lty = 1,
  txt.adj = NULL,
  txt.pos = 2,
  txt.cex = 0.8,
  txt.col = "black"
)

Arguments

Details

The function allows users to vary certain plotting parameters to produce different graphical outcomes for plotRefToTarget. Not all parameters need to be adjusted to use this function, as the defaults above will be used.

Author(s)

Michael Collyer & Emma Sherratt

See Also

plotRefToTarget, text, text3d

Examples

## Not run: 
data(plethodon) 
Y.gpa <- gpagen(plethodon$land)    #GPA-alignment
ref <- mshape(Y.gpa$coords)
plotRefToTarget(ref, Y.gpa$coords[,,39]) # default settings

# Altering points and links
GP1 <- gridPar(pt.bg = "red", pt.size = 1, link.col = "blue", link.lwd = 2, 
n.col.cell = 50)
plotRefToTarget(ref, Y.gpa$coords[,,39], gridPars = GP1, mag = 2, 
links = plethodon$links, method = "TPS")

# Altering point color
GP2 <- gridPar(pt.bg = "green", pt.size = 1) 
plotRefToTarget(ref, Y.gpa$coords[,,39], gridPars = GP2, mag = 3, 
method = "vector")

# Altering ref and target points
GP3 <- gridPar(pt.bg = "blue", pt.size = 1.5, tar.pt.bg = "orange", 
tar.pt.size = 1) 
plotRefToTarget(ref, Y.gpa$coords[,,39], gridPars = GP3, mag = 3, 
method = "points")

# Altering outline color
GP4 <- gridPar(tar.out.col = "red", tar.out.cex = 0.3) 
plotRefToTarget(ref, Y.gpa$coords[,,39], gridPars = GP4, mag = 3, 
outline = plethodon$outline, method = "TPS")

# Altering text labels
GP5 <- gridPar(txt.pos = 3, txt.col = "red") 
plotRefToTarget(ref, Y.gpa$coords[,,39], gridPars = GP5, mag = 3, 
method = "vector", label = TRUE)

## End(Not run)

Description

Landmark data from hummingbird bills (includes sliding semilandmarks on curves)

Author(s)

Chelsea Berns and Dean Adams

References

Berns, C.M., and Adams, D.C. 2010. Bill shape and sexual shape dimorphism between two species of temperate hummingbirds: Archilochus alexandri (black-chinned hummingbirds) and Archilochus colubris (ruby-throated hummingbirds). The Auk. 127:626-635.

Description

Function quantifies the degree of morphological integration between modules of Procrustes shape variables

Usage

integration.test(
  A,
  A2 = NULL,
  partition.gp = NULL,
  iter = 999,
  seed = NULL,
  print.progress = TRUE
)

Arguments

Details

The function quantifies the degree of morphological integration between modular partitions of shape data as defined by Procrustes shape variables. It is assumed that the landmarks have previously been aligned using Generalized Procrustes Analysis (GPA) [e.g., with gpagen]. The function may be used to assess the degree of morphological integration between two or more sets of variables.

The function estimates the degree of morphological integration using a two-block partial least squares analysis (PLS). When used with landmark data, this analysis is referred to as singular warps analysis (Bookstein et al. 2003). If more than two partitions are defined, the average pairwise PLS correlation is utilized as the test statistic. The observed test value is then compared to a distribution of values obtained by randomly permuting the individuals (rows) in one partition relative to those in the other. A significant result is found when the observed PLS correlation is large relative to this distribution, and implies that the structures are integrated with one another (see Bookstein et al. 2003). In addition, a multivariate effect size describing the strength of the effect is estimated from the empirically-generated sampling distribution (see details in Adams and Collyer 2019). If only two partitions are specified, a plot of PLS scores along the first set of PLS axes is optionally displayed, and thin-plate spline deformation grids along these axes are also shown if data were input as a 3D array.

Input for the analysis can take one of two forms. First, one can input a single dataset (as a matrix or 3D array, along with a vector describing which variables correspond to which partitions (for the case of a 3D array, which landmarks belong to which partitions is specified). Alternatively, when evaluating the integration between two structures or partitions, two datasets may be provided.

The generic functions, print, summary, and plot all work with integration.test. The generic function, plot, produces a two-block.pls plot. This function calls plot.pls, which produces an ordination plot. An additional argument allows one to include a vector to label points. Starting with version 3.1.0, warpgrids are no longer available with plot.pls but after making a plot, the function returns values that can be used with picknplot.shape or a combination of shape.predictor and plotRefToTarget to visualize shape changes in the plot (via warpgrids).

Similarity to two.b.pls and compare.pls

Note that integration.test performed on two matrices or arrays returns the same results as two.b.pls. However, two.b.pls is limited to only two modules. It might be of interest with 3+ modules to perform integration tests between all pairwise comparisons of modules. This can be done, test by test, and the levels of integration can be compared with compare.pls. Such results are different than using the average amount of integration, as performed by integration.test when more than two modules are input.

Using phylogenies and PGLS

If one wishes to incorporate a phylogeny, phylo.integration is the function to use. This function is exactly the same as integration.test but allows PGLS estimation of PLS vectors. Because integration.test can be used on two blocks, phylo.integration likewise allows one to perform a phylogenetic two-block PLS analysis.

Notes for geomorph 3.0.4 and subsequent versions

Compared to previous versions of geomorph, users might notice differences in effect sizes. Previous versions used z-scores calculated with expected values of statistics from null hypotheses (sensu Collyer et al. 2015); however Adams and Collyer (2016) showed that expected values for some statistics can vary with sample size and variable number, and recommended finding the expected value, empirically, as the mean from the set of random outcomes. Geomorph 3.0.4 and subsequent versions now center z-scores on their empirically estimated expected values and where appropriate, log-transform values to assure statistics are normally distributed. This can result in negative effect sizes, when statistics are smaller than expected compared to the average random outcome. For ANOVA-based functions, the option to choose among different statistics to measure effect size is now a function argument.

Value

Objects of class "pls" from integration.test return a list of the following:

Author(s)

Dean Adams

References

Bookstein, F. L., P. Gunz, P. Mitteroecker, H. Prossinger, K. Schaefer, and H. Seidler. 2003. Cranial integration in Homo: singular warps analysis of the midsagittal plane in ontogeny and evolution. J. Hum. Evol. 44:167-187.

Collyer, M.L., D.J. Sekora, and D.C. Adams. 2015. A method for analysis of phenotypic change for phenotypes described by high-dimensional data. Heredity. 115:357-365.

Adams, D.C. and M.L. Collyer. 2016. On the comparison of the strength of morphological integration across morphometric datasets. Evolution. 70:2623-2631.

Adams, D.C. and M.L. Collyer. 2019. Comparing the strength of modular signal, and evaluating alternative modular hypotheses, using covariance ratio effect sizes with morphometric data. Evolution. 73:2352-2367.

See Also

two.b.pls, modularity.test, phylo.integration, and compare.pls

Examples

## Not run: 
data(plethodon) 
Y.gpa <- gpagen(plethodon$land)    #GPA-alignment   
 
#landmarks on the skull and mandible assigned to partitions
land.gps <- c("A","A","A","A","A","B","B","B","B","B","B","B") 

IT <- integration.test(Y.gpa$coords, partition.gp = land.gps)
summary(IT) # Test summary
P <- plot(IT) # PLS plot

 # Visualize shape at minimum and maximum PLS scores.
 # This is the challenging way
 
 # Block 1 
 minx <- min(P$plot_args$x)
 maxx <- max(P$plot_args$x)
 preds <- shape.predictor(P$A1, 
 x = P$plot.args$x,
 min = minx, max = maxx)
 plotRefToTarget(mshape(P$A1), preds$min)
 plotRefToTarget(mshape(P$A1), preds$max)

 # Block 2 
 miny <- min(P$plot_args$y)
 maxy <- max(P$plot_args$y)
 preds <- shape.predictor(P$A2, 
 x = P$plot.args$y,
 min = miny, max = maxy)
 plotRefToTarget(mshape(P$A2), preds$min)
 plotRefToTarget(mshape(P$A2), preds$max)
 
 ### Visualize shape variation using picknplot.shape Because picknplot  
 ### requires user decisions, the following example
 ### is not run (but can be with removal of #).
 ### For detailed options, see the picknplot help file
 # picknplot.shape(P)

IT$left.pls.vectors # extracting just the left (first block) 
# singular vectors

## End(Not run)

Description

Function quantifies the morphological integration in a set of traits

Usage

integration.Vrel(A, phy = NULL)

Arguments

Details

The function quantifies the strength of morphological integration in a set of variables. Here the set of traits are treated as a single unit, and the overall degree of covariation in them is quantified using the relative eigenvalue index: Vrel (Pavlicev et al. 2009). Following Conaway and Adams (2022), only the non-trivial dimensions of variation are used in the calculation of Vrel. The measure is then converted to an effect size (Z-score), based on the procedures in Conaway and Adams (2022). These may be used in subsequent comparisons of the strength of integration across datasets. Input for the analysis may be a 3D array of Procrustes coordinates, of a matrix of variables. If the observations are species related by a phylogeny, the phylogeny may also be included.

Value

Objects of class "rel.eig" from integration.Vrel return a list of the following:

Author(s)

Dean Adams

References

Pavlicev, M., J. M. Cheverud, and G. P. Wagner. 2009. Measuring morphological integration using eigenvalue variance. Evolutionary Biology 36:157-170.

Conaway, M.A., and D.C. Adams. 2022. An effect size for comparing the strength of morphological integration across studies. Evolution. 76: 2244-2259.

See Also

compare.ZVrel

Examples

## Not run: 
data(plethodon) 
Y.gpa <- gpagen(plethodon$land)    #GPA-alignment    
integration.Vrel(Y.gpa$coords)

## End(Not run)

Description

A function to calculate linear distances between a set of landmark coordinates (interlandmark distances)

Usage

interlmkdist(A, lmks)

Arguments

Details

Function takes a 3D array of landmark coordinates from a set of specimens and the addresses for the start and end landmarks defining linear measurements and then calculates the interlandmark distances. The function returns a matrix of linear distances for all specimens. If the 'lmks' matrix has row or column names defining the name of the linear measurements, the returned matrix will use these for column names (see example). If only two interlandmark distances, 'lmks' input must be m x 2.

Value

Function returns a matrix (n x m) of m linear distances for n specimens

Author(s)

Emma Sherratt & Michael Collyer

Examples

## Not run: 
data(plethodon)
# Make a matrix defining three interlandmark distances 
lmks <- matrix(c(8,9,6,12,4,2), ncol=2, byrow=TRUE, 
dimnames = list(c("eyeW", "headL", "mouthL"),c("start", "end")))

# where 8-9 is eye width; 6-12 is head length; 4-2 is mouth length
# or alternatively

lmks <- data.frame(eyeW = c(8,9), headL = c(6,12), mouthL = c(4,2), 
row.names = c("start", "end")) 

A <- plethodon$land
lineardists <- interlmkdist(A, lmks)

## End(Not run)

Description

Landmark data from heads and tails of larval Salamanders exposed to different treatments of herbicides.

Details

Data set includes tail landmarks (coords), index for identifying semilandmarks (sliders), and vectors for variables including herbicide treatment (Treatment) and Family (Family). The latter variable indicates the egg masses (clutches) sampled in the wild from which individual eggs were randomly assigned to treatment. See Levis et al. (2016) for more experimental details.

Author(s)

Michael Collyer

References

Levis, N.A, M.L. Schooler, J.R. Johnson, and M.L. Collyer. 2016. The effects of terrestrial and aquatic herbicides on larval salamander morphology and swim speed. Biological Journal of the Linnean Society. 118:569-581.

Description

Superimposed landmark data of the dorsal view of lizard heads

Details

Dataset includes superimposed landmarks (coords), centroid size (cs), an index of individuals (ind) and digitizing repetitions (rep), and a table of symmetrical matching landmarks (lm.pairs). The object is a geomorph.data.frame. The dataset corresponds to the data for population "b" from Lazic et al. 2015.

Author(s)

Antigoni Kaliontzopoulou

References

Lazic, M., Carretero, M.A., Crnobrnja-Isailovic, J. & Kaliontzopoulou, A. 2015. Effects of environmental disturbance on phenotypic variation: an integrated assessment of canalization, developmental stability, modularity and allometry in lizard head shape. American Naturalist 185: 44-58.

Description

Function attempts to coerce plot information from a geomorph plot object to an amenable ggplot object.

Usage

make_ggplot(object)

Arguments

Details

This function will attempt to use the plot arguments from an geomorph plot object to make a ggplot that can be additionally updated, as desired. Not all plot characteristics might be converted. Nonetheless, a ggplot will be coerced and could be updated, according to user preference.

This function assumes no responsibility for arguments made by ggplot. It merely produces a ggplot object that should resemble a geomorph plot default. Any augmentation of ggplot objects can be done either by direct intervention of the ggplot produced or reformatting the initial geomorph plot produced. One should not expect direct correspondence between R base plot parameters and ggplot parameters.

Author(s)

Michael Collyer

Examples

## Not run: 

### PLS Example
 data(plethodon) 
 Y.gpa <- gpagen(plethodon$land)    #GPA-alignment    
 landmarks on the skull and mandible assigned to partitions
 land.gps <- c("A","A","A","A","A","B","B","B","B","B","B","B") 
 IT <- integration.test(Y.gpa$coords, partition.gp = land.gps, iter = 999)
 summary(IT) # Test summary
 P <- plot(IT) # PLS plot
 make_ggplot(P) # same plot in ggplot

### Allometry example

 data(plethodon) 
 Y.gpa <- gpagen(plethodon$land, print.progress = FALSE)    #GPA-alignment  

 gdf <- geomorph.data.frame(Y.gpa, site = plethodon$site, 
                           species = plethodon$species) 

 fit <- procD.lm(coords ~ Csize * species * site, data=gdf, iter=0, 
                 print.progress = FALSE)

 P <- plotAllometry(fit, size = gdf$Csize, logsz = TRUE, method = "PredLine", 
                     pch = 19, col = as.numeric(interaction(gdf$species, gdf$site)))

 make_ggplot(P)

### Tangent Space plot

 data(plethspecies) 
 Y.gpa <- gpagen(plethspecies$land)    #GPA-alignment

 PCA.w.phylo <- gm.prcomp(Y.gpa$coords, phy = plethspecies$phy)
 P <- plot(PCA.w.phylo, phylo = TRUE, main = "PCA.w.phylo")
 make_ggplot(P)

## End(Not run)

Description

Function quantifies the degree of modularity between two or more hypothesized modules of Procrustes shape variables and compares this to patterns found by randomly assigning landmarks into subsets

Usage

modularity.test(
  A,
  partition.gp,
  iter = 999,
  CI = FALSE,
  seed = NULL,
  opt.rot = TRUE,
  print.progress = TRUE
)

Arguments

Details

The function quantifies the degree of modularity in two or more hypothesized modules of Procrustes shape variables, and compares this to what is expected under the null hypothesis of random assignment of variables to partitions (i.e., neither modular nor integrated structure). Input may be either a 2D matrix of phenotypic values, or a 3D array of Procrustes shape variables. It is assumed that the landmarks have previously been aligned using Generalized Procrustes Analysis (GPA); e.g., with gpagen. The degree of modularity is quantified using the CR coefficient (Adams 2016). If more than two modules are defined, the average pairwise CR coefficient is utilized. The CR coefficient for the observed modular hypothesis is then compared to a distribution of values obtained by randomly assigning landmarks into subsets, with the restriction that the number of landmarks in each subset is identical to that observed in each of the original partitions. A significant modular signal is found when the observed CR coefficient is small relative to this distribution (see Adams 2016). Such a result implies that there is significantly greater independence among modules than is expected under the null hypothesis of random associations of variables (neither modular nor integrated structure). This result is consistent with the identification of significant modular structure in the data. In addition, a multivariate effect size describing the strength of the effect is estimated from the empirically-generated sampling distribution (see details in Adams and Collyer 2019). A histogram of coefficients obtained via resampling is presented, with the observed value designated by an arrow in the plot. For landmark data, the CR coefficient found from the average CR across a 90 degree rotation of the data is used as the test statistic (see Adams 2016). For all data, the CR coefficient is returned, and (optionally) its 95 for landmark data, estimation of the CI can take some time to compute.

Landmark groups can be defined using define.modules, or made by hand (see example below). To use this method with other data (i.e., a set of length measurements), the input A should be a matrix of n rows of specimens and variables arranged in columns. In this case, the partition.gp input should have each variable assigned to a partition.

The generic functions, print, summary, and plot all work with modularity.test. The generic function, plot, produces a histogram of random CR values associated with the resampling procedure.

Value

An object of class "CR" is a list containing the following

Author(s)

Dean Adams

References

Adams, D.C. 2016.Evaluating modularity in morphometric data: Challenges with the RV coefficient and a new test measure. Methods in Ecology and Evolution 7:565-572.

Adams, D.C. and M.L. Collyer. 2019. Comparing the strength of modular signal, and evaluating alternative modular hypotheses, using covariance ratio effect sizes with morphometric data. Evolution. 73:2352-2367.

See Also

two.b.pls, integration.test, phylo.modularity, and phylo.integration

Examples

## Not run: 
 data(pupfish) 
 Y.gpa <- gpagen(pupfish$coords, print.progress = FALSE)    #GPA-alignment    
 #landmarks on the body and operculum
 land.gps <- rep('a',56); land.gps[39:48] <- 'b'

 MT <- modularity.test(Y.gpa$coords, land.gps, CI = FALSE)
 summary(MT) # Test summary
 plot(MT) # Histogram of CR sampling distribution 
# Result implies modularity present

## End(Not run)

Description

Function estimates morphological disparity and performs pairwise comparisons among groups

Usage

morphol.disparity(
  f1,
  groups = NULL,
  partial = FALSE,
  transform = FALSE,
  iter = 999,
  seed = NULL,
  data = NULL,
  print.progress = TRUE,
  ...
)

Arguments

Details

The function estimates morphological disparity and performs pairwise comparisons to identify differences among groups. Morphological disparity is estimated as the Procrustes variance, overall or for groups, using residuals of a linear model fit. Procrustes variance is the same sum of the diagonal elements of the group covariance matrix divided by the number of observations in the group (e.g., Zelditch et al. 2012). The function takes as input a formula to describe the linear model fit, plus a formulaic indication of groups (e.g., ~ groups). It is assumed that the formula describes shape data that have been GPA-aligned [e.g., gpagen], although the function can work with any multivariate data.

Partial disparities (Foote 1993) can also be calculated, but only for model formulas containing only an intercept (e.g., coords ~ 1). Partial disparity has the same numerator as Procrustes variance (with respect to an overall mean) but the denominator is N - 1 for all N observations, rather than n, the group size. (The sum of all group n equals N.) Partial disparities have the appeal that the sum of group partial disparities it the total disparity.

Absolute differences in Procrustes variances are test statistics that can be used to test differences in morphological disparity among groups. These differences are statistically evaluated through permutation, where the vectors of residuals are randomized among groups. The function can be used to obtain disparity for the whole dataset by using "a dummy group factor "~ 1" as the right-hand portion of the formula, in which case only Procrustes variance is returned. Additionally, if the right-hand portion of the formula only contains (continuous) covariates, e.g., "~ Csize", Procrustes variance will be calculated for the whole data set or groups, after accounting for the linear regression described. Finally, different factors can be indicated in the formula and for groups, if one wishes to compare morphological disparities for groups comprising only a portion of or collapsing of the groups in a more complex model (see examples).

This function can be used with an object of class "procD.lm" or "lm.rrpp", if such analyses have already been performed. This is especially useful for analyses performed with procD.pgls. In this case, residuals obtained from PGLS estimation of coefficients, rather than OLS estimation, will be used in the analysis. Thus, one can account for phylogeny when comparing morphological disparity among groups. However, one should be aware that PGLS finds the PGLS residuals and transforms them via the phylogenetic covariance matrix in order to calculate variance. Thus, disparity (dispersion in tangent space) and variance (dispersion in transformed space, after accounting for phylogeny) can be considered different things. This function uses an argument, transform, to allow users to use either GLS residuals in tangent space or transformed residuals in the transformed space to calculate disparity. The former characterizes the dispersion of actual shapes around GLS means; the latter estimates GLS variances, by group.

Notes for geomorph 3.1.0 and subsequent versions

The function pairwise in the RRPP package can also be used to evaluate morphological disparity, and will yield results identical to those of the current function. A simple example is shown below

Value

Objects of class "morphol.disparity" return a list with the following components (if groups are specified):

Author(s)

Emma Sherratt and Michael Collyer

References

Zelditch, M. L., D. L. Swiderski, H. D. Sheets, and W. L. Fink. 2012. Geometric morphometrics for biologists: a primer. 2nd edition. Elsevier/Academic Press, Amsterdam.

Foote, M. 1993. Contributions of individual taxa to overall morphological disparity. Paleobiology, 19: 403-419.

Collyer, M.L and D.C. Adams, 2021. Phylogenetically Aligned component analysis. Methods in Ecology and Evolution, 12: 369-372.

Examples

## Not run: 

 data(plethodon)
 Y.gpa <- gpagen(plethodon$land, print.progress = FALSE)    #GPA-alignment
 gdf <- geomorph.data.frame(Y.gpa, species = plethodon$species, 
 site = plethodon$site)

# Morphological disparity for entire data set
 morphol.disparity(coords ~ 1, groups = NULL, data = gdf, 
 print.progress = FALSE)

# Morphological disparity for entire data set, accounting for allometry
 morphol.disparity(coords ~ Csize, groups= NULL, data = gdf, 
 print.progress = FALSE)

# Morphological disparity without covariates, using overall mean
 morphol.disparity(coords ~ 1, groups= ~ species * site, data = gdf, 
 print.progress = FALSE)

# Morphological partial disparities for overal mean
 morphol.disparity(coords ~ 1, groups= ~ species * site, partial = TRUE, 
 data = gdf, print.progress = FALSE)

# Morphological disparity without covariates, using group means
 morphol.disparity(coords ~ species * site, groups= ~species * site, 
 data = gdf, print.progress = FALSE)

# Morphological disparity of different groups than those 
# described by the linear model
 morphol.disparity(coords ~ Csize + species * site, groups= ~ species, 
 data = gdf, print.progress = FALSE)

# Extracting components
 MD <- morphol.disparity(coords ~ Csize + species * site, groups= ~ species, 
 data = gdf, print.progress = FALSE)
 MD$Procrustes.var # just the Procrustes variances


### Morphol.disparity can be used with previously-defined 
### procD.lm or lm.rrpp class objects

 data(plethspecies)
 Y.gpa <- gpagen(plethspecies$land)    #GPA-alignment
 gp.end <- factor(c(0,0,1,0,0,1,1,0,0))  #endangered species vs. rest
 names(gp.end) <- plethspecies$phy$tip

 gdf <- geomorph.data.frame(Y.gpa, phy = plethspecies$phy, 
 gp.end = gp.end)

 pleth.ols <- procD.lm(coords ~ Csize + gp.end, 
 data = gdf) # ordinary least squares
 pleth.pgls <- procD.pgls(coords ~ Csize + gp.end, phy = phy, 
 data = gdf) # phylogenetic generalized least squares

 summary(pleth.ols)
 summary(pleth.pgls)

 morphol.disparity(f1 = pleth.ols, groups = ~ gp.end, data = gdf, 
 print.progress = FALSE)
 morphol.disparity(f1 = pleth.pgls, groups = ~ gp.end, 
 transform = FALSE, data = gdf, 
 print.progress = FALSE) # disparity in tangent space
 morphol.disparity(f1 = pleth.pgls, groups = ~ gp.end,
 transform = TRUE, data = gdf, 
 print.progress = FALSE) # disparity in transformed space 

# Three plots that correspond to the three tests
 PCA <- gm.prcomp(Y.gpa$coords, phy = plethspecies$phy)
 pPCA <- gm.prcomp(Y.gpa$coords, phy = plethspecies$phy, 
 GLS = TRUE, transform = FALSE)
 tpPCA <- gm.prcomp(Y.gpa$coords, phy = plethspecies$phy, 
 GLS = TRUE, transform = TRUE)

 par(mfrow = c(1,3))

# Phylomorphospace
 PC.plot <- plot(PCA, pch = 19, phylo = TRUE, main = "PCA-OLS")
 shapeHulls(PC.plot, groups = gp.end)

# Phylo-PCA
 pPC.plot <- plot(pPCA, pch = 19, phylo = TRUE, main = "pPCA - GLS, not transformed")
 shapeHulls(pPC.plot, groups = gp.end)

# Transformed phylo-PCA
 tpPC.plot <- plot(tpPCA, pch = 19, phylo = TRUE, main = "tpPCA - GLS, transformed")
 shapeHulls(tpPC.plot, groups = gp.end)

 par(mfrow = c(1,1))

### Variance using RRPP (not necessarily the same as disparity)
 PW <- pairwise(pleth.ols, groups = gp.end)
 summary(PW, test.type = 'var')
 PW2 <- pairwise(pleth.pgls, groups = gp.end)
 summary(PW2, test.type = 'var')

## End(Not run)

Description

Landmarks on mosquito wings

Author(s)

Dean Adams

Description

Estimate the mean shape for a set of aligned specimens

Usage

mshape(A, na.action = 1)

Arguments

Details

The function estimates the average landmark coordinates for a set of aligned specimens. Three different methods are available for missing data (see Arguments and Examples).

One can then use the generic function plot to produce a numbered plot of landmark positions and potentially add links, in order to review landmark positions

Author(s)

Antigoni Kaliontzopoulou & Michael Collyer

Examples

## Not run: 
data(plethodon)
Y.gpa <- gpagen(plethodon$land)    #GPA-alignment
A <- Y.gpa$coords
A[[1]] <- NA # make a missing value, just for example

mshape(Y.gpa$coords)   # mean (consensus) configuration
# mshape(A, na.action = 1) # will return an error
mshape(A, na.action = 2) # returns NA in spot of missing value
mshape(A, na.action = 3) # finds mean values from all possible values

## End(Not run)

Description

Handle missing values in rrpp.data.frame objects

Usage

## S3 method for class 'geomorph.data.frame'
na.omit(object, ...)

Arguments

Author(s)

Michael Collyer

Examples

data(plethspecies)
Y.gpa <- gpagen(plethspecies$land, verbose = TRUE)
gdf <- geomorph.data.frame(Y.gpa)
gdf$d <- Y.gpa$procD
gdf$group <- c(rep(1, 4), rep(2, 4), NA) # one unknown group designation
gdf
ndf <- na.omit(gdf)
ndf

Description

Function quantifies the degree of phylogenetic morphological covariation between two or more sets of Procrustes shape variables using partial least squares.

Usage

phylo.integration(
  A,
  A2 = NULL,
  phy,
  partition.gp = NULL,
  iter = 999,
  seed = NULL,
  print.progress = TRUE
)

Arguments

Details

The function quantifies the degree of phylogenetic morphological integration between two or more sets of Procrustes shape variables. The approach is based on a Brownian motion model of evolution. It is assumed that the landmarks have previously been aligned using Generalized Procrustes Analysis (GPA) [e.g., with gpagen].

The function estimates the degree of morphological covariation between two or sets of variables while accounting for phylogeny using partial least squares (Adams and Felice 2014), and under a Brownian motion model of evolution. If more than two partitions are defined, the average pairwise PLS correlation is utilized as the test statistic. The observed value is statistically assessed using permutation, where data for one partition are permuted relative to the other partitions. In addition, a multivariate effect size describing the strength of the effect is estimated from the empirically-generated sampling distribution (see details in Adams and Collyer 2019). Note that this permutation is performed on phylogenetically- transformed data, so that the probability of phylogenetic association of A vs. B is similar to that of B vs. A: i.e., prob(A,B|phy)~prob(B,A|phy); thus, shuffling the correct exchangeable units under the null hypothesis of no integration (Adams and Collyer 2018).

Input for the analysis can take one of two forms. First, one can input a single dataset (as a matrix or 3D array, along with a vector describing which variables correspond to which partitions (for the case of a 3D array, which landmarks belong to which partitions is specified). Alternatively, when evaluating the integration between two structures or partitions, two datasets may be provided.

The generic functions, print, summary, and plot all work with phylo.integration. The generic function, plot, produces a two-block.pls plot. This function calls plot.pls, which produces an ordination plot. An additional argument allows one to include a vector to label points. Starting with version 3.1.0, warpgrids are no longer available with plot.pls but after making a plot, the function returns values that can be used with picknplot.shape or a combination of shape.predictor and plotRefToTarget to visualize shape changes in the plot (via warpgrids).

Similarity to two.b.pls and compare.pls

Note that phylo.integration performed on two matrices or arrays returns the same results as a phylogenetic variation of two.b.pls. It might be of interest with 3+ modules to perform separate phylogenetic integration tests between all pairwise comparisons of modules. This can be done, test by test, and the levels of integration can be compared with compare.pls. Such results are different than using the average amount of integration when more than two modules are input, as found with phylo.integration.

Notes for geomorph 3.0.4 and subsequent versions

Compared to previous versions of geomorph, users might notice differences in effect sizes. Previous versions used z-scores calculated with expected values of statistics from null hypotheses (sensu Collyer et al. 2015); however Adams and Collyer (2016) showed that expected values for some statistics can vary with sample size and variable number, and recommended finding the expected value, empirically, as the mean from the set of random outcomes. Geomorph 3.0.4 and subsequent versions now center z-scores on their empirically estimated expected values and where appropriate, log-transform values to assure statistics are normally distributed. This can result in negative effect sizes, when statistics are smaller than expected compared to the average random outcome. For ANOVA-based functions, the option to choose among different statistics to measure effect size is now a function argument.

Value

Objects of class "pls" from integration.test return a list of the following:

Author(s)

Dean Adams

References

Adams, D.C. and R. Felice. 2014. Assessing phylogenetic morphological integration and trait covariation in morphometric data using evolutionary covariance matrices. PLOS ONE. 9(4):e94335.

Collyer, M.L., D.J. Sekora, and D.C. Adams. 2015. A method for analysis of phenotypic change for phenotypes described by high-dimensional data. Heredity. 115:357-365.

Adams, D.C. and M.L. Collyer. 2016. On the comparison of the strength of morphological integration across morphometric datasets. Evolution. 70:2623-2631.

Adams, D.C. and M.L. Collyer. 2018. Multivariate comparative methods: evaluations, comparisons, and recommendations. Systematic Biology. 67:14-31.

Adams, D.C. and M.L. Collyer. 2019. Comparing the strength of modular signal, and evaluating alternative modular hypotheses, using covariance ratio effect sizes with morphometric data. Evolution. 73:2352-2367.

See Also

integration.test, modularity.test, and two.b.pls

Examples

## Not run: 
data(plethspecies) 
Y.gpa <- gpagen(plethspecies$land)    #GPA-alignment
land.gps <- c("A","A","A","A","A","B","B","B","B","B","B") 

IT <- phylo.integration(Y.gpa$coords, partition.gp = land.gps,
  phy = plethspecies$phy)
summary(IT) # Test summary
P <- plot(IT) # PLS plot

 # Block 1 
 minx <- min(P$plot_args$x)
 maxx <- max(P$plot_args$x)
 preds <- shape.predictor(P$A1, 
 x = P$plot.args$x,
 min = minx, max = maxx)
 plotRefToTarget(mshape(P$A1), preds$min)
 plotRefToTarget(mshape(P$A1), preds$max)

 # Block 2 
 miny <- min(P$plot_args$y)
 maxy <- max(P$plot_args$y)
 preds <- shape.predictor(P$A2, 
 x = P$plot.args$y,
 min = miny, max = maxy)
 plotRefToTarget(mshape(P$A2), preds$min)
 plotRefToTarget(mshape(P$A2), preds$max)
 
 ### Visualize shape variation using picknplot.shape Because picknplot  
 ### requires user decisions, the following example
 ### is not run.
 ### For detailed options, see the picknplot help file
 # picknplot.shape(P)

## End(Not run)

Description

Function quantifies the degree of modularity between two or more hypothesized modules of Procrustes shape variables in a phylogenetic context and compares this to patterns found by randomly assigning landmarks into subsets

Usage

phylo.modularity(
  A,
  partition.gp,
  phy,
  CI = FALSE,
  iter = 999,
  seed = NULL,
  print.progress = TRUE
)

Arguments

Details

The function quantifies the degree of phylogenetic modularity in two or more hypothesized modules of Procrustes shape variables as defined by landmark coordinates, under a Brownian motion model of evolution. The degree of modularity is characterized by the covariance ratio covariance ratio (CR: see Adams 2016). The phylogenetic version of the approach procedure utilizes the evolutionary covariance matrix among traits found under a Brownian motion model of evolution as the basis of the analysis. This is the same matrix used to evaluate patterns of phylogenetic morphological integration as described in Adams and Felice (2014).

Input may be either a 2D matrix of phenotypic values, or a 3D array of Procrustes shape variables. It is assumed that the landmarks have previously been aligned using Generalized Procrustes Analysis (GPA) [e.g., with gpagen]. The degree of modularity is quantified using the CR coefficient (Adams 2016). If more than two modules are defined, the average pairwise CR coefficient is utilized. The CR coefficient for the observed modular hypothesis is then compared to a distribution of values obtained by randomly assigning landmarks into subsets, with the restriction that the number of landmarks in each subset is identical to that observed in each of the original partitions. A significant modular signal is found when the observed CR coefficient is small relative to this distribution (see Adams 2016). Such a result implies that there is significantly greater independence among modules than is expected under the null hypothesis of random associations of variables (neither modular nor integrated structure). This result is consistent with the identification of significant modular structure in the data. For landmark data, the CR coefficient found from the average CR across a 90 degree rotation of the data is used as the test statistic (see Adams 2016). In addition, a multivariate effect size describing the strength of the effect is estimated from the empirically-generated sampling distribution (see details in Adams and Collyer 2019).

Value

Objects of class "CR" from modularity.test return a list of the following:

Author(s)

Dean Adams

References

Adams, D.C. 2016.Evaluating modularity in morphometric data: Challenges with the RV coefficient and a new test measure. Methods in Ecology and Evolution. 7:565-572.

Adams, D.C. and R. Felice. 2014. Assessing phylogenetic morphological integration and trait covariation in morphometric data using evolutionary covariance matrices. PLOS ONE. 9(4):e94335.

Adams, D.C. and M.L. Collyer. 2019. Comparing the strength of modular signal, and evaluating alternative modular hypotheses, using covariance ratio effect sizes with morphometric data. Evolution. 73:2352-2367.

Examples

## Not run: 

 data(plethspecies)
 Y.gpa <- gpagen(plethspecies$land)    #GPA-alignment
 land.gps <- c("A","A","A","A","A","B","B","B","B","B","B") 

 MT <- phylo.modularity(Y.gpa$coords, partition.gp = land.gps, 
 phy = plethspecies$phy, 
 CI = FALSE)
 summary(MT) # Test summary
 plot(MT) # Histogram of CR sampling distribution 
 
## End(Not run)

Description

Function calculates the degree of phylogenetic signal from Procrustes shape variables

Usage

physignal(A, phy, iter = 999, seed = NULL, print.progress = FALSE)

Arguments

Details

The function estimates the degree of phylogenetic signal present in Procrustes shape variables for a given phylogeny. It is assumed that the landmarks have previously been aligned using Generalized Procrustes Analysis (GPA) [e.g., with gpagen]. The degree of phylogenetic signal in data is estimated using the multivariate version of the K-statistic (Kmult: Adams 2014). This value evaluates the degree of phylogenetic signal in a dataset relative to what is expected under a Brownian motion model of evolution. For geometric morphometric data, the approach is a mathematical generalization of the Kappa statistic (Blomberg et al. 2003) appropriate for highly multivariate data (see Adams 2014). Significance testing is found by permuting the shape data among the tips of the phylogeny. And effect size can be estimated based on likelihood method that finds and optimal scaling of the tree; see physignal.z.

This function can also be used with univariate data (i.e. centroid size) if imported as matrix with rownames giving the taxa names. In this case, the estimate of phylogenetic signal is identical to that found using the standard kappa statistic (Blomberg et al. 2003).

The generic functions, print, summary, and plot all work with physignal. The generic function, plot, produces a histogram of random K statistics, associated with the resampling procedure.

Notes for geomorph 4.0

Starting with version 4.0.2, this function no longer reports an effect size. Effect sizes for phylogenetic signal can be calculated with physignal.z, based on likelihood.

Notes for geomorph 3.0

Compared to older versions of geomorph, the order of input variables has changed, so that it is consistent with other functions in the program. Additionally, users should note that the function physignal no longer contains multiple methods. Only Kmult is used. Thus, for older scripts method="" should be removed from the function call.

Value

Function returns a list with the following components:

Author(s)

Dean Adams and Michael Collyer

References

Blomberg SP, Garland T, Ives AR. 2003. Testing for phylogenetic signal in comparative data: behavioral traits are more labile. Evolution, 57:717-745.

Adams, D.C. 2014. A generalized K statistic for estimating phylogenetic signal from shape and other high-dimensional multivariate data. Systematic Biology. 63:685-697.

Adams, D.C. and M.L. Collyer. 2019. Comparing the strength of modular signal, and evaluating alternative modular hypotheses, using covariance ratio effect sizes with morphometric data. Evolution. 73:2352-2367.

See Also

gm.prcomp, physignal.z

Examples

## Not run: 
data(plethspecies) 
Y.gpa <- gpagen(plethspecies$land)    #GPA-alignment    

#Test for phylogenetic signal in shape
PS.shape <- physignal(A = Y.gpa$coords, phy = plethspecies$phy)
summary(PS.shape)
plot(PS.shape)
plot(PS.shape$PACA, phylo = TRUE)
PS.shape$K.by.p # Phylogenetic signal profile

#Test for phylogenetic signal in size
PS.size <- physignal(A = Y.gpa$Csize, phy = plethspecies$phy)
summary(PS.size)
plot(PS.size)

## End(Not run)

Description

Function calculates the degree of phylogenetic signal across dimensions of shape space from Procrustes shape variables

Usage

physignal.eigen(
  Y,
  phy = NULL,
  Cov = NULL,
  Blomberg = FALSE,
  unit.tree = TRUE,
  lambda = 1,
  test = TRUE,
  iter = 999,
  seed = NULL,
  tol = 0.001
)

Arguments

Details

The function estimates the degree of phylogenetic signal across dimensions of shape space, based on Procrustes shape variables for a given phylogeny. The method is intended to interrogate patterns of phylogenetic signal that are not uniformly distributed across shape space, but rather are concentrated in one or a few dimensions. The degree of phylogenetic signal in data is estimated using a matrix generalization of Blomberg's Kappa (K), and the phylogenetic signal represented by this matrix is decomposed into linear combinations that describe directions of phylogenetic signal within the dataspace containing decreasing levels of phylogenetic signal (Mitteroecker et al. 2024). Two summary measures based on the eigenvalues of K are calculated: the determinant of K (detK) and the trace of K (traceK). Both describe the degree of phylogenetic signal in multivariate data. In addition, the multivariate version of the K-statistic (Kmult: Adams 2014) is also provided. All three statistics are evaluated via permutation, where the shape data is permuted among the tips of the phylogeny. Effect sizes are also provided.

It should be noted that by default the function uses OLS centering when calculating the phylogenetic signal matrix, unlike most authors who have used GLS centering. Both OLS and GLS centering were proposed in Blomberg et al.s (2003) original exploration of phylogenetic signal. However, for both mathematical and computational reasons OLS is used here (see justification in Mitteroecker et al. 2024). Both measures are highly rank correlated with one another (Mitteroecker et al. 2024). Additionally, using using BLOMBERG = TRUE will result in GLS centering, in which case the statistic Kmult obtained with this function will be identical to that obtained when using physignal.

Importantly, because detK and traceK are based on the covariance matrix K, singularity can become an issue. In particular, geometric morphometric shape data will not be of full rank, as several dimensions are standardized during the Generalized Procrustes Analysis (one may also have fewer observations than variables, which will also generate redundancies). For this reason, a principal components analysis of the data is performed, and the redundant dimensions are removed so that detK and traceK may be computed (see Mitteroecker et al. 2024). Additionally, if n< (p X k), the last nontrivial PC dimension is also removed, as in this case, using 100 variation results in invariant K-statistics across permutations.

The generic functions, print, summary, and plot all work with physignal.eigen.

Value

Function returns a list with the following components:

Author(s)

Dean Adams and Michael Collyer

References

Mitteroecker, P., M.L. Collyer, and D.C. Adams. 2024. Exploring phylogenetic signal in multivariate phenotypes by maximizing Blomberg's K. Systematic Biology. (In Press).

Blomberg SP, Garland T, Ives AR. 2003. Testing for phylogenetic signal in comparative data: behavioral traits are more labile. Evolution, 57:717-745.

Adams, D.C. 2014. A generalized K statistic for estimating phylogenetic signal from shape and other high-dimensional multivariate data. Systematic Biology. 63:685-697.

See Also

gm.prcomp, physignal

Examples

## Not run: 
data(plethspecies) 
Y.gpa <- gpagen(plethspecies$land)    #GPA-alignment    

#Test for phylogenetic signal in shape
PSe.shape <- physignal.eigen(Y = Y.gpa$coords, phy = plethspecies$phy)
summary(PSe.shape)
plot(PSe.shape)
plot(PSe.shape, type = "vectors")


## End(Not run)

Description

Function calculates the degree of phylogenetic signal from Procrustes shape variables, using effect size estimated from likelihood.

Usage

physignal.z(
  A,
  phy,
  lambda = c("burn", "mean", "front", "all"),
  iter = 999,
  seed = NULL,
  tol = 1e-04,
  PAC.no = NULL,
  print.progress = FALSE,
  verbose = FALSE
)