Addition by Subtraction: Improving geomorph with fewer functions

Michael Collyer

2023-01-16

Introduction

The evolution of an R package includes both additions and deletions. Some functions remain necessary, some are replaced by better functions, and others prove either problematic or redundant if there are better alternatives. Over the last few years, geomorph has mostly evolved by additions. The repertoire of functions is now rather impressive. However, recently, some of the functions have become redundant, especially as functions in the package, RRPP, have offered better alternatives for the analytical functions in geomorph (specifically, those that use randomization of residuals in permutation procedures, RRPP). Therefore, with version 3.1.0, we decided to deprecate some functions in lieu of better options available in RRPP. This vignette illustrates how to use alternative functions and update scripts for veterans, but should also be a good introduction for new users.

The following is a summary of major changes in version 3.1.0, showing old functions and new alternatives:

geomorph 3.0.7 function status geomorph 3.1.0 / RRPP alternative
advanced.procD.lm deprecated procD.lm + pairwise (RRPP)
procD.allometry deprecated procD.lm + anova (RRPP) + plotAllometry (new geomorph function)
nested.update deprecated anova (RRPP) (Can do much more than nested.update could)
trajectory.analysis moved Now in RRPP with improved functionality

By eliminating or updating some of the older functions, the overall functionality of geomorph has increased (addition by subtraction). To illustrate the improvements in detail, we focus on one data example below.

Data example: larval Salamanders

These data consist of head landmarks, tail landmarks, and indicators for their semilandmarks (sliders) for larval salamanders, from an experimental study (Levis et al. 2016). Salamanders from different families (egg clutches) were exposed to different herbicide treatments. For our examples, we will focus on tail shape. First, let’s perform GPA.

library(geomorph)
data("larvalMorph")

Y.gpa <- gpagen(larvalMorph$tailcoords, curves = larvalMorph$tail.sliders,
                ProcD = FALSE, print.progress = FALSE)
plot(Y.gpa)

Next, let’s define some putative models of shape variation that would be biologically relevant.

gdf <- geomorph.data.frame(Y.gpa, treatment = larvalMorph$treatment, 
                           family = larvalMorph$family)

fit.size <- procD.lm(coords ~ log(Csize), 
                     data = gdf, print.progress = FALSE) # simple allometry model
fit.family<- procD.lm(coords ~ log(Csize) * family, 
                     data = gdf, print.progress = FALSE) # unique family allometries
fit.treatment<- procD.lm(coords ~ log(Csize) * treatment/family, 
                     data = gdf, print.progress = FALSE) # unique treatment: family allometries

fit.size
## 
## Linear Model fit with lm.rrpp
## 
## Number of observations: 114
## Number of dependent variables: 128  
## Data space dimensions: 113  
## Sums of Squares and Cross-products: Type I
## Number of permutations: 1000
## Call: procD.lm(f1 = coords ~ log(Csize), data = gdf, print.progress = FALSE)
fit.family
## 
## Linear Model fit with lm.rrpp
## 
## Number of observations: 114
## Number of dependent variables: 128  
## Data space dimensions: 113  
## Sums of Squares and Cross-products: Type I
## Number of permutations: 1000
## Call: procD.lm(f1 = coords ~ log(Csize) * family, data = gdf, print.progress = FALSE)
fit.treatment
## 
## Linear Model fit with lm.rrpp
## 
## Number of observations: 114
## Number of dependent variables: 128  
## Data space dimensions: 113  
## Sums of Squares and Cross-products: Type I
## Number of permutations: 1000
## Call: procD.lm(f1 = coords ~ log(Csize) * treatment/family, data = gdf,  
##     print.progress = FALSE)

In geomorph 3.0.7, one could use summary to produce an ANOVA table. That still works in geomorph 3.1.0, but the anova function (in RRPP) is much more flexible!

How to perform ANOVA and model comparisons in geomorph 3.1.0

Replaces advanced.procD.lm and enhances functionality

The following are the typical ANOVA tables for each model:

anova(fit.size)
## 
## Analysis of Variance, using Residual Randomization
## Permutation procedure: Randomization of null model residuals 
## Number of permutations: 1000 
## Estimation method: Ordinary Least Squares 
## Sums of Squares and Cross-products: Type I 
## Effect sizes (Z) based on F distributions
## 
##             Df      SS      MS     Rsq      F      Z Pr(>F)   
## log(Csize)   1 0.50263 0.50263 0.16206 21.661 3.1141  0.001 **
## Residuals  112 2.59885 0.02320 0.83794                        
## Total      113 3.10147                                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Call: procD.lm(f1 = coords ~ log(Csize), data = gdf, print.progress = FALSE)
anova(fit.family)
## 
## Analysis of Variance, using Residual Randomization
## Permutation procedure: Randomization of null model residuals 
## Number of permutations: 1000 
## Estimation method: Ordinary Least Squares 
## Sums of Squares and Cross-products: Type I 
## Effect sizes (Z) based on F distributions
## 
##                    Df      SS      MS     Rsq       F       Z Pr(>F)   
## log(Csize)          1 0.50263 0.50263 0.16206 22.5073  3.1333  0.001 **
## family              3 0.21009 0.07003 0.06774  3.1360  2.2172  0.011 * 
## log(Csize):family   3 0.02160 0.00720 0.00696  0.3224 -1.5886  0.952   
## Residuals         106 2.36715 0.02233 0.76324                          
## Total             113 3.10147                                          
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Call: procD.lm(f1 = coords ~ log(Csize) * family, data = gdf, print.progress = FALSE)
anova(fit.treatment)
## 
## Analysis of Variance, using Residual Randomization
## Permutation procedure: Randomization of null model residuals 
## Number of permutations: 1000 
## Estimation method: Ordinary Least Squares 
## Sums of Squares and Cross-products: Type I 
## Effect sizes (Z) based on F distributions
## 
##                              Df      SS      MS     Rsq       F       Z Pr(>F)
## log(Csize)                    1 0.50263 0.50263 0.16206 26.6673  3.2716  0.001
## treatment                     5 0.45807 0.09161 0.14770  4.8607  3.5710  0.001
## log(Csize):treatment          5 0.08025 0.01605 0.02587  0.8515 -0.1177  0.553
## log(Csize):treatment:family  18 0.47729 0.02652 0.15389  1.4068  1.3236  0.098
## Residuals                    84 1.58323 0.01885 0.51048                       
## Total                       113 3.10147                                       
##                               
## log(Csize)                  **
## treatment                   **
## log(Csize):treatment          
## log(Csize):treatment:family . 
## Residuals                     
## Total                         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Call: procD.lm(f1 = coords ~ log(Csize) * treatment/family, data = gdf,  
##     print.progress = FALSE)

In geomorph 3.0.7, one could use advanced.procD.lm to perform ANOVA on one reduced and one full model. With anova in RRPP, one can directly compare any number of models. For example

anova(fit.size, fit.family, fit.treatment, print.progress = FALSE)
## 
## Analysis of Variance, using Residual Randomization
## Permutation procedure: Randomization of null model residuals 
## Number of permutations: 1000 
## Estimation method: Ordinary Least Squares 
## Effect sizes (Z) based on F distributions
## 
##                                        ResDf Df    RSS      SS       MS     Rsq
## coords ~ log(Csize) (Null)               112  1 2.5989                  0.00000
## coords ~ log(Csize) * family             106  6 2.3672 0.23169 0.038616 0.07470
## coords ~ log(Csize) * treatment/family    84 28 1.5832 1.01562 0.036272 0.32746
## Total                                    113    3.1015                         
##                                             F      Z     P Pr(>F)
## coords ~ log(Csize) (Null)                                       
## coords ~ log(Csize) * family           1.7292 1.4262 0.083       
## coords ~ log(Csize) * treatment/family 1.9244 2.7124 0.005       
## Total

One might ascertain that models that include parameters for family and treatment are significant improvements over a model that only contains size. Like advanced.procD.lm, one can perform model comparisons, but now one can do that on any number of models with one model serving as the null model for comparisons. Note that the model fits can use types I, II, or III, sums of squares and cross-products (SSCP), ordinary or generalized least squares, and the ANOVA can be performed on a number of different test statistics (see anova.lm.rrpp help file for more details).

How to perform allometry analyses

Alternatives to procD.allometry and advanced.procD.lm, with increased flexibility

The following analyses replace procD.allometry, primarily, and reiterate the enhanced method of multiple model comparisons. First, let’s just consider the simple allometry model

summary(fit.size)
## 
## Analysis of Variance, using Residual Randomization
## Permutation procedure: Randomization of null model residuals 
## Number of permutations: 1000 
## Estimation method: Ordinary Least Squares 
## Sums of Squares and Cross-products: Type I 
## Effect sizes (Z) based on F distributions
## 
##             Df      SS      MS     Rsq      F      Z Pr(>F)   
## log(Csize)   1 0.50263 0.50263 0.16206 21.661 3.1141  0.001 **
## Residuals  112 2.59885 0.02320 0.83794                        
## Total      113 3.10147                                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Call: procD.lm(f1 = coords ~ log(Csize), data = gdf, print.progress = FALSE)

It is clear that there is a significant association between shape and size. To visualize this relationship, there are many options. First, we can use the plot generic for procD.lm and vary some of the arguments. Let’s do that, using two different ways to visualize shape change: prediction lines (PredLine) and regression scores (RegScore). The former are first principal component scores for fitted values from the procD.lm fit (Adams and Nistri 2010); the latter are standardized projected shape scores, along the axis defined by the regression of shape on size (Drake and Klingenberg 2008). These plot options were formerly found in plot.procD.allometry, as well as plot.procD.lm (as shown below).

plot(fit.size, type = "regression", reg.type = "PredLine", predictor = log(gdf$Csize))

plot(fit.size, type = "regression", reg.type = "RegScore", predictor = log(gdf$Csize))

The plot.procD.lm function is for any procD.lm fit. If one wishes to work specifically with allometry models, the plotAllometry function performs the same analysis as a convenient wrapper for plot.procD.lm. For example,

plotAllometry(fit.size, size = gdf$Csize, logsz = TRUE, method = "PredLine")

plotAllometry(fit.size, size = gdf$Csize, logsz = TRUE, method = "RegScore")

An important detail with these plots is that PredLine and RegScore are model-based projections of shape data. As we will see below, changing the model changes the outcome of the plot.

We could also perform a two-block partial least squares (PLS) analysis to find the correlation between shape and size, which is not based on a particular model.

PLS <- two.b.pls(log(gdf$Csize), gdf$coords, print.progress = FALSE)
PLS
## 
## Call:
## two.b.pls(A1 = log(gdf$Csize), A2 = gdf$coords, print.progress = FALSE) 
## 
## 
## 
## r-PLS: 0.466
## 
## Effect Size (Z): 3.1676
## 
## P-value: 0.001
## 
## Based on 1000 random permutations
plot(PLS)

An astute observer might catch that the PLS plot is exactly the same as the RegScore plot. It is, in this case of a simple allometry model. They are also both the same as a plot of the common allometric component (CAC, Mitteroecker et al. 2004); i.e.,

plotAllometry(fit.size, size = gdf$Csize, logsz = TRUE, method = "CAC")

The CAC plot will always be the same as the PLS plot, irrespective of the type of shape-allometry model. The RegScore plot is the same in this simple case because only one vector of regression coefficients is produced, which aligns perfectly with the major axis of covariation between shape and size (the CAC or the shape PLS vector; Adams et al. 2013)

One can also append a size vector to a matrix of shape variables and perform principal components analysis (PCA), called size-shape PCA (Mitteroecker et al. 2004).

plotAllometry(fit.size, size = gdf$Csize, logsz = TRUE, method = "size.shape")

More complex allometry models

We already learned that family and treatment were “significant” model effects. Now let’s focus on whether we should believe that families or treatments have unique allometries or a common allometry

fit.unique <- procD.lm(coords ~ log(Csize) * treatment/family, 
                     data = gdf, print.progress = FALSE) # unique allometries
fit.common <- procD.lm(coords ~ log(Csize) + treatment/family, 
                     data = gdf, print.progress = FALSE) # common allometry
anova(fit.common, fit.unique, print.progress = FALSE)
## 
## Analysis of Variance, using Residual Randomization
## Permutation procedure: Randomization of null model residuals 
## Number of permutations: 1000 
## Estimation method: Ordinary Least Squares 
## Effect sizes (Z) based on F distributions
## 
##                                               ResDf Df    RSS       SS       MS
## coords ~ log(Csize) + treatment/family (Null)    89  1 1.6587                  
## coords ~ log(Csize) * treatment/family           84  5 1.5832 0.075471 0.015094
## Total                                           113    3.1015                  
##                                                    Rsq      F        Z     P
## coords ~ log(Csize) + treatment/family (Null) 0.000000                      
## coords ~ log(Csize) * treatment/family        0.024334 0.8008 -0.23384 0.598
## Total                                                                       
##                                               Pr(>F)
## coords ~ log(Csize) + treatment/family (Null)       
## coords ~ log(Csize) * treatment/family              
## Total

The ANOVA above was formerly carried out by advanced.procD.lm as a homogeneity of slopes (HOS) test in the function procD.allometry (although nested effects were not possible.) Because this model comparison did not yield a significant result, we can conclude that a common allometry model is appropriate. Thus, we might want to plot the results, color-coding the points by treatment

plotAllometry(fit.common, size = gdf$Csize, logsz = TRUE, method = "PredLine",
              pch = 19, col = as.numeric(gdf$treatment))

plotAllometry(fit.common, size = gdf$Csize, logsz = TRUE, method = "RegScore",
              pch = 19, col = as.numeric(gdf$treatment))

The next section focuses on ANOVA for model effects and pairwise comparisons.

How to update ANOVA and perform pairwise comparisons in geomorph 3.1.0

Replaces advanced.procD.lm and nested.update, and enhances functionality

In our example, we have both fixed and random effects. Treatment is a fixed effect and family is a random effect, nested within our fixed effect (as egg clutches were “randomly”” sampled from the wild). Generally, when evaluating model effects, ANOVA involves assessing the probability of observed F-values which are ratios of mean squared (MS) values for effects to MS values for the appropriate random effect, usually the residuals. (For procD.lm models, the distribution of F-values is generated over many random permutations.) For example,

anova(fit.common)
## 
## Analysis of Variance, using Residual Randomization
## Permutation procedure: Randomization of null model residuals 
## Number of permutations: 1000 
## Estimation method: Ordinary Least Squares 
## Sums of Squares and Cross-products: Type I 
## Effect sizes (Z) based on F distributions
## 
##                   Df      SS      MS     Rsq       F      Z Pr(>F)   
## log(Csize)         1 0.50263 0.50263 0.16206 26.9691 3.2807  0.001 **
## treatment          5 0.45807 0.09161 0.14770  4.9157 3.5979  0.001 **
## treatment:family  18 0.48207 0.02678 0.15543  1.4370 1.4376  0.075 . 
## Residuals         89 1.65870 0.01864 0.53481                         
## Total            113 3.10147                                         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Call: procD.lm(f1 = coords ~ log(Csize) + treatment/family, data = gdf,  
##     print.progress = FALSE)

Notice that the F-value for the three effects - log(Csize), Treatment, and Treatment:Family - is calculated as MS effect / MS Residuals. This is the default. However, in our mixed-model ANOVA, we would prefer to calculate the F-value for treatment as MS Treatment / MS Treatment:Family, to determine if the treatment effect is meaningful compared to shape variation among families we sampled randomly and assigned to treatments. We can update our ANOVA by specifying what the denominator (error term) should be, as a sequence of error terms for the ANOVA; e.g.,

anova(fit.common, error = c("Residuals", "treatment:family", "Residuals"))
## 
## Analysis of Variance, using Residual Randomization
## Permutation procedure: Randomization of null model residuals 
## Number of permutations: 1000 
## Estimation method: Ordinary Least Squares 
## Sums of Squares and Cross-products: Type I 
## Effect sizes (Z) based on F distributions
## 
##                   Df      SS      MS     Rsq       F      Z Pr(>F)   
## log(Csize)         1 0.50263 0.50263 0.16206 26.9691 3.2807  0.001 **
## treatment          5 0.45807 0.09161 0.14770  3.4208 1.9482  0.025 * 
## treatment:family  18 0.48207 0.02678 0.15543  1.4370 1.4376  0.075 . 
## Residuals         89 1.65870 0.01864 0.53481                         
## Total            113 3.10147                                         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Call: procD.lm(f1 = coords ~ log(Csize) + treatment/family, data = gdf,  
##     print.progress = FALSE)

Notice the F-value and effect size decreased a bit with recalculation of the F-value for treatment, but remained significant. What this function did was recalculate every treatment F-value in every random permutation to generate a new distribution for measuring effect size (Z) and P-value.

This same procedure was accomplished with the nested.update function in geomorph 3.0.7, but was limited to one term nested within another. The anova function with error argument presented here can work on any number of random effects and interactions, allowing much more complex models to be evaluated!

Pairwise comparisons to determine which treatments differ

Now that we know that shape covaries with size, but in a common way for each treatment, we might wish to compare treatment least-squares (LS) means to see which treatments differ in shape, accounting for allometry and accounting for family effects. In geomorph 3.0.7, this was accomplished with advanced.procD.lm. In geomorph 3.1.0, we use the pairwise function from RRPP. There are several advantages to using the pairwise function. The most prominent are: (1) the ability to quickly change among different test types and summaries without re-analysis, (2) alternative summary options, (3) an option to use grouping variables not included in the original model fit, and (4) much faster computation for straightforward tests.

The pairwise function has this general format:

pairwise(fit, groups, covariate),

where fit is an already fitted object, using procD.lm, groups is a factor to designate groups to be compared, and covariate is a vector if slopes are to be compared. This format assumes the inherent null model of “fit” is appropriate. If an alternative null model is desired, the function can be updated as:

pairwise(fit, fit.null, groups, covariate),

where fit.null is a second procD.lm fit. If one is not sure about the inherent null model, they can use the reveal.model.designs function of RRPP to discover the exact null model used; e.g.,

reveal.model.designs(fit.common)
##                                  Reduced 
## log(Csize)                             1 
## treatment                     log(Csize) 
## treatment:family  log(Csize) + treatment 
##                                                       Full
## log(Csize)                                      log(Csize)
## treatment                           log(Csize) + treatment
## treatment:family log(Csize) + treatment + treatment:family
##                                                   
## log(Csize)                                        
## treatment                                         
## treatment:family <- Null/Full inherent in pairwise

The results tell us that if we run pairwise on fit.common, the null model would be ~ log(Csize) + treatment and the full model would be ~ log(Csize) + treatment + treatment:family. This is the case because we used type I (sequential) sums of squares and cross-products (the default). However, it is maybe not ideal. We might prefer to have as a null model, ~ log(Csize) + family. Thus, let’s first establish that model and then run the pairwise function

fit.null <- procD.lm(coords ~ log(Csize) + family, data = gdf, print.progress = FALSE)
PW <- pairwise(fit.common, fit.null, groups = gdf$treatment, print.progress = FALSE)
PW
## 
## Pairwise comparisons
## 
## Groups: AG AM CF CI TG TM 
## 
## RRPP: 1000 permutations

There are now many options for summarizing results; i.e., we can perform multiple tests! Here is one option:

summary(PW, test.type = "dist", confidence = 0.95)
## 
## Pairwise comparisons
## 
## Groups: AG AM CF CI TG TM 
## 
## RRPP: 1000 permutations
## 
## LS means:
## Vectors hidden (use show.vectors = TRUE to view)
## 
## Pairwise distances between means, plus statistics
##                d  UCL (95%)          Z Pr > d
## AG:AM 0.03031748 0.08544041 -0.3983502  0.625
## AG:CF 0.06775829 0.08787927  1.1355437  0.144
## AG:CI 0.12685901 0.09140712  2.3133960  0.004
## AG:TG 0.10993275 0.08320799  2.1315610  0.010
## AG:TM 0.06447743 0.08823741  1.0581303  0.162
## AM:CF 0.05742092 0.08919549  0.8048851  0.227
## AM:CI 0.12904332 0.08977658  2.3407942  0.006
## AM:TG 0.11418302 0.08878564  2.1302365  0.010
## AM:TM 0.04829385 0.08787553  0.5319077  0.309
## CF:CI 0.13108751 0.08447949  2.5457646  0.003
## CF:TG 0.13207622 0.08781678  2.4905214  0.003
## CF:TM 0.06377460 0.09022502  0.9788791  0.190
## CI:TG 0.05102561 0.09150320  0.5093922  0.324
## CI:TM 0.17096322 0.09526279  2.7628577  0.001
## TG:TM 0.15667496 0.08594155  2.8008853  0.001

The test statistics used, “dist”, is the distance between LS means. By specifying a confidence level, we are given upper confidence limits (UCL) from the distributions of pairwise distances. We can see that if the observed distance is larger than the UCL, the P-value is less than 1 - confidence; i.e., it is “significant”. The default is this “stats table”, but we could also produce pairwise tables. In fact, we can reproduce the old format for advanced.procD.lm like so:

anova(fit.null, fit.common, print.progress = FALSE)
## 
## Analysis of Variance, using Residual Randomization
## Permutation procedure: Randomization of null model residuals 
## Number of permutations: 1000 
## Estimation method: Ordinary Least Squares 
## Effect sizes (Z) based on F distributions
## 
##                                        ResDf Df    RSS      SS       MS     Rsq
## coords ~ log(Csize) + family (Null)      109  1 2.3887                  0.00000
## coords ~ log(Csize) + treatment/family    89 20 1.6587 0.73005 0.036503 0.23539
## Total                                    113    3.1015                         
##                                             F     Z     P Pr(>F)
## coords ~ log(Csize) + family (Null)                             
## coords ~ log(Csize) + treatment/family 1.9586 2.463 0.011       
## Total
summary(PW, test.type = "dist", confidence = 0.95, stat.table = FALSE)
## 
## Pairwise comparisons
## 
## Groups: AG AM CF CI TG TM 
## 
## RRPP: 1000 permutations
## 
## LS means:
## Vectors hidden (use show.vectors = TRUE to view)
## 
## Pairwise distances between means
##            AG         AM         CF         CI         TG         TM
## AG 0.00000000 0.03031748 0.06775829 0.12685901 0.10993275 0.06447743
## AM 0.03031748 0.00000000 0.05742092 0.12904332 0.11418302 0.04829385
## CF 0.06775829 0.05742092 0.00000000 0.13108751 0.13207622 0.06377460
## CI 0.12685901 0.12904332 0.13108751 0.00000000 0.05102561 0.17096322
## TG 0.10993275 0.11418302 0.13207622 0.05102561 0.00000000 0.15667496
## TM 0.06447743 0.04829385 0.06377460 0.17096322 0.15667496 0.00000000
## 
## Pairwise 95% Upper confidence limits between means
##            AG         AM         CF         CI         TG         TM
## AG 0.00000000 0.08544041 0.08787927 0.09140712 0.08320799 0.08823741
## AM 0.08544041 0.00000000 0.08919549 0.08977658 0.08878564 0.08787553
## CF 0.08787927 0.08919549 0.00000000 0.08447949 0.08781678 0.09022502
## CI 0.09140712 0.08977658 0.08447949 0.00000000 0.09150320 0.09526279
## TG 0.08320799 0.08878564 0.08781678 0.09150320 0.00000000 0.08594155
## TM 0.08823741 0.08787553 0.09022502 0.09526279 0.08594155 0.00000000
## 
## Pairwise effect sizes (Z) between means
##            AG         AM        CF        CI        TG        TM
## AG  0.0000000 -0.3983502 1.1355437 2.3133960 2.1315610 1.0581303
## AM -0.3983502  0.0000000 0.8048851 2.3407942 2.1302365 0.5319077
## CF  1.1355437  0.8048851 0.0000000 2.5457646 2.4905214 0.9788791
## CI  2.3133960  2.3407942 2.5457646 0.0000000 0.5093922 2.7628577
## TG  2.1315610  2.1302365 2.4905214 0.5093922 0.0000000 2.8008853
## TM  1.0581303  0.5319077 0.9788791 2.7628577 2.8008853 0.0000000
## 
## Pairwise P-values between means
##       AG    AM    CF    CI    TG    TM
## AG 1.000 0.625 0.144 0.004 0.010 0.162
## AM 0.625 1.000 0.227 0.006 0.010 0.309
## CF 0.144 0.227 1.000 0.003 0.003 0.190
## CI 0.004 0.006 0.003 1.000 0.324 0.001
## TG 0.010 0.010 0.003 0.324 1.000 0.001
## TM 0.162 0.309 0.190 0.001 0.001 1.000

Because we have already performed the pairwise procedure, we could also summarize a different test. For example, let’s say we wish to compare morphological disparities (variances) among treatments. We simply change the summary:

summary(PW, test.type = "var", confidence = 0.95)
## 
## Pairwise comparisons
## 
## Groups: AG AM CF CI TG TM 
## 
## RRPP: 1000 permutations
## 
## 
## Observed variances by group
## 
##         AG         AM         CF         CI         TG         TM 
## 0.01419375 0.01384366 0.01923172 0.01704117 0.01006673 0.01287420 
## 
## Pairwise distances between variances, plus statistics
##                  d   UCL (95%)           Z Pr > d
## AG:AM 0.0003500900 0.010087309 -1.69635391  0.949
## AG:CF 0.0050379659 0.009938407  0.52168706  0.305
## AG:CI 0.0028474188 0.010372592 -0.23118837  0.590
## AG:TG 0.0041270191 0.010316444  0.22999925  0.433
## AG:TM 0.0013195570 0.009850599 -0.87460939  0.789
## AM:CF 0.0053880560 0.009613089  0.59377947  0.300
## AM:CI 0.0031975089 0.010284072 -0.06700465  0.548
## AM:TG 0.0037769290 0.009932049  0.19034887  0.443
## AM:TM 0.0009694669 0.009679284 -1.08080047  0.850
## CF:CI 0.0021905471 0.010327023 -0.48165241  0.674
## CF:TG 0.0091649850 0.010154643  1.40254488  0.075
## CF:TM 0.0063575229 0.010012406  0.85125709  0.211
## CI:TG 0.0069744379 0.010679392  0.89506258  0.202
## CI:TM 0.0041669758 0.010341725  0.26148124  0.410
## TG:TM 0.0028074621 0.010485985 -0.20477972  0.578

This should be exactly the same as performing a morphological disparity test

morphol.disparity(fit.common, groups = gdf$treatment, print.progress = FALSE)
## 
## Call:
## morphol.disparity(f1 = fit.common, groups = gdf$treatment, print.progress = FALSE) 
## 
## 
## 
## 
## Randomized Residual Permutation Procedure Used
## 1000 Permutations
## 
## Procrustes variances for defined groups
##         AG         AM         CF         CI         TG         TM 
## 0.01419375 0.01384366 0.01923172 0.01704117 0.01006673 0.01287420 
## 
## 
## Pairwise absolute differences between variances
##             AG           AM          CF          CI          TG           TM
## AG 0.000000000 0.0003500900 0.005037966 0.002847419 0.004127019 0.0013195570
## AM 0.000350090 0.0000000000 0.005388056 0.003197509 0.003776929 0.0009694669
## CF 0.005037966 0.0053880560 0.000000000 0.002190547 0.009164985 0.0063575229
## CI 0.002847419 0.0031975089 0.002190547 0.000000000 0.006974438 0.0041669758
## TG 0.004127019 0.0037769290 0.009164985 0.006974438 0.000000000 0.0028074621
## TM 0.001319557 0.0009694669 0.006357523 0.004166976 0.002807462 0.0000000000
## 
## 
## P-Values
##       AG    AM    CF    CI    TG    TM
## AG 1.000 0.949 0.305 0.590 0.433 0.789
## AM 0.949 1.000 0.300 0.548 0.443 0.850
## CF 0.305 0.300 1.000 0.674 0.075 0.211
## CI 0.590 0.548 0.674 1.000 0.202 0.410
## TG 0.433 0.443 0.075 0.202 1.000 0.578
## TM 0.789 0.850 0.211 0.410 0.578 1.000

The pairwise function in RRPP is really versatile and far less cumbersome than advanced.procD.lm was. More examples are provided in the help file for the function, and greater detail for how to summarize different tests is found in the summary.pairwise help file.

How to perform trajectory analysis in geomorph 3.1.0

Uses trajectory.analysis, now in RRPP but formerly in geomorph

The trajectory,analysis function has been modified and now resides in RRPP, as it is a wrapper for the pairwise function. It has the same basic arguments as the pairwise function, but also has an argument for trajectory points (which can be a single value, if data are already trajectories or a factor to indicate trajectory point levels). Following the example above, trajectory analysis can be considered a pairwise function where treatments are trajectories and families are trajectory points. The following highlights the steps involved for one type of example (but the plotting options are quite numerous):

TA <- trajectory.analysis(fit.common, 
                          groups = gdf$treatment, traj.pts = gdf$family,
                          pca = TRUE, print.progress = FALSE)
summary(TA, attribute = "MD")
## 
## Trajectory analysis
## 
## 1000 permutations.
## 
## Points projected onto trajectory PCs
## 
## Trajectories:
## Trajectories hidden (use show.trajectories = TRUE to view)
## 
## Observed path distances by group
## 
##        AG        AM        CF        CI        TG        TM 
## 0.2122330 0.1916723 0.6845214 0.2411477 0.2356911 0.1719081 
## 
## Pairwise absolute differences in path distances, plus statistics
##                d UCL (95%)          Z Pr > d
## AG:AM 0.02056069 0.2260201 -0.9724758  0.825
## AG:CF 0.47228842 0.2315721  2.9165995  0.001
## AG:CI 0.02891473 0.2373102 -0.7759162  0.780
## AG:TG 0.02345814 0.2335531 -0.9983408  0.828
## AG:TM 0.04032491 0.2185481 -0.5317521  0.701
## AM:CF 0.49284911 0.2341886  2.9900560  0.001
## AM:CI 0.04947542 0.2364577 -0.3521631  0.633
## AM:TG 0.04401883 0.2308622 -0.4905885  0.684
## AM:TM 0.01976423 0.2392499 -1.0665493  0.842
## CF:CI 0.44337369 0.2338955  2.6931482  0.001
## CF:TG 0.44883028 0.2484048  2.7824283  0.001
## CF:TM 0.51261333 0.2368642  3.1037237  0.001
## CI:TG 0.00545659 0.2304833 -1.8682112  0.969
## CI:TM 0.06923965 0.2477956 -0.1031573  0.536
## TG:TM 0.06378306 0.2456721 -0.0810123  0.537
TP <- plot(TA, pch = 19, cex = 0.7, col = as.numeric(gdf$treatment))
add.trajectories(TP, traj.bg = 1:nlevels(gdf$treatment), 
                 start.bg = 1:nlevels(gdf$treatment),
                 end.bg = 1:nlevels(gdf$treatment))

The argument, attribute = “MD”, indicates that the differences between trajectory magnitudes - the path length of trajectories - is considered. The trajectory analysis could also be summarized for trajectory correlations (TC), which are angles between trajectory directions (major axes of variation), or shape differences (SD). More examples are given in the trajectory.analysis help file. The plotting options for trajectory.analysis in geomorph 3.0.7 were somewhat constrained but have been enhanced in RRPP. The function, plot.trajectory.analysis plots the data points projected on the PCs for fitted values and the function, add.trajectories, superimposes the trajectories on these points. The help files for these functions have further details.

Summary

This vignette hopefully illustrates that the goals of former functions - advanced.procD.lm, procD.allometry, and nested.update - are as easily accomplished with alternative functions that have fewer constraints. If you are new to geomorph, then this vignette should give you some ideas for how to analyze data with pairwise comparisons.

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.

Adams, D. C., and A. Nistri. 2010. Ontogenetic convergence and evolution of foot morphology in European cave salamanders (Family: Plethodontidae). BMC Evol. Biol. 10:1-10.

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.

Drake, A. G., and C. P. Klingenberg. 2008. The pace of morphological change: Historical transformation of skull shape in St Bernard dogs. Proc. R. Soc. B. 275:71-76.

Mitteroecker, P., P. Gunz, M. Bernhard, K. Schaefer, and F. L. Bookstein. 2004. Comparison of cranial ontogenetic trajectories among great apes and humans. J. Hum. Evol. 46:679-698.

Rohlf, F.J., and M. Corti. 2000. The use of partial least-squares to study covariation in shape. Systematic Biology 49: 740-753.