Addition by Subtraction: Improving geomorph with fewer functions

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	deprectaed	anova (RRPP) (Can do much more than nested.update could)
trajectory.analysis	moved	Now in RRPP with improved functionality

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)

## Loading required package: RRPP

## Loading required package: rgl

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.43502 0.43502 0.1623 21.699 3.605  0.001 **
## Residuals  112 2.24536 0.02005 0.8377                       
## Total      113 2.68038                                      
## ---
## 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.43502 0.43502 0.16230 22.5517  3.6377  0.001 **
## family              3 0.18019 0.06006 0.06722  3.1136  2.1784  0.010 * 
## log(Csize):family   3 0.02044 0.00681 0.00763  0.3532 -1.5400  0.943   
## Residuals         106 2.04473 0.01929 0.76285                          
## Total             113 2.68038                                          
## ---
## 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
## log(Csize)                    1 0.43502 0.43502 0.16230 26.9257  3.8040
## treatment                     5 0.41046 0.08209 0.15314  5.0811  3.5517
## log(Csize):treatment          5 0.07121 0.01424 0.02657  0.8815 -0.0507
## log(Csize):treatment:family  18 0.40655 0.02259 0.15168  1.3980  1.3310
## Residuals                    84 1.35713 0.01616 0.50632                
## Total                       113 2.68038                                
##                             Pr(>F)   
## log(Csize)                   0.001 **
## treatment                    0.001 **
## log(Csize):treatment         0.548   
## log(Csize):treatment:family  0.099 . 
## 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
## coords ~ log(Csize) (Null)               112  1 2.2454                 
## coords ~ log(Csize) * family             106  6 2.0447 0.20063 0.033438
## coords ~ log(Csize) * treatment/family    84 28 1.3571 0.88822 0.031722
## Total                                    113    2.6804                 
##                                            Rsq      F      Z     P Pr(>F)
## coords ~ log(Csize) (Null)             0.16230                           
## coords ~ log(Csize) * family           0.23715 1.7334 1.8854 0.022       
## coords ~ log(Csize) * treatment/family 0.49368 1.9635 3.8061 0.001       
## Total

One might ascertain that models that include paramaters 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 SS distributions
## 
##             Df      SS      MS    Rsq      F      Z Pr(>SS)   
## log(Csize)   1 0.43502 0.43502 0.1623 21.699 3.4672   0.001 **
## Residuals  112 2.24536 0.02005 0.8377                         
## Total      113 2.68038                                        
## ---
## 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 correaltion 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.475
## 
## 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
## coords ~ log(Csize) + treatment/family (Null)    89  1 1.4221         
## coords ~ log(Csize) * treatment/family           84  5 1.3571 0.064998
## Total                                           113    2.6804         
##                                                  MS     Rsq      F      Z
## coords ~ log(Csize) + treatment/family (Null)       0.46943              
## coords ~ log(Csize) * treatment/family        0.013 0.49368 0.8046 1.0723
## Total                                                                    
##                                                   P Pr(>F)
## coords ~ log(Csize) + treatment/family (Null)             
## coords ~ log(Csize) * treatment/family        0.134       
## 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 conlude 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.43502 0.43502 0.16230 27.2245 3.8203  0.001 **
## treatment          5 0.41046 0.08209 0.15314  5.1375 3.5884  0.001 **
## treatment:family  18 0.41276 0.02293 0.15399  1.4351 1.4601  0.070 . 
## Residuals         89 1.42213 0.01598 0.53057                         
## Total            113 2.68038                                         
## ---
## 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 Treament: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.43502 0.43502 0.16230 27.2245 3.8203  0.001 **
## treatment          5 0.41046 0.08209 0.15314  3.5799 2.0989  0.016 * 
## treatment:family  18 0.41276 0.02293 0.15399  1.4351 1.4601  0.070 . 
## Residuals         89 1.42213 0.01598 0.53057                         
## Total            113 2.68038                                         
## ---
## 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 eaxt 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.03003071 0.07675692 -0.4761783  0.594
## AG:CF 0.06598651 0.08148254  1.2593319  0.121
## AG:CI 0.11537143 0.08379440  3.5878927  0.005
## AG:TG 0.10498004 0.07699149  3.4659402  0.006
## AG:TM 0.06321247 0.08157740  1.1375433  0.130
## AM:CF 0.05611306 0.08216896  0.7376608  0.202
## AM:CI 0.11643446 0.08273533  3.5504394  0.007
## AM:TG 0.10810704 0.08113110  3.4251391  0.005
## AM:TM 0.04798806 0.08009783  0.3995548  0.262
## CF:CI 0.12400279 0.07738504  4.4056219  0.003
## CF:TG 0.13013212 0.07915375  4.6313544  0.001
## CF:TM 0.06024489 0.08217325  0.9170473  0.183
## CI:TG 0.04780469 0.08263148  0.2730329  0.325
## CI:TM 0.15838941 0.08728829  5.1622430  0.001
## TG:TM 0.15012798 0.07862142  5.6782533  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
## coords ~ log(Csize) + family (Null)      109  1 2.0652                 
## coords ~ log(Csize) + treatment/family    89 20 1.4221 0.64304 0.032152
## Total                                    113    2.6804                 
##                                            Rsq      F      Z     P Pr(>F)
## coords ~ log(Csize) + family (Null)    0.22952                           
## coords ~ log(Csize) + treatment/family 0.46943 2.0121 3.8573 0.001       
## 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.03003071 0.06598651 0.11537143 0.10498004 0.06321247
## AM 0.03003071 0.00000000 0.05611306 0.11643446 0.10810704 0.04798806
## CF 0.06598651 0.05611306 0.00000000 0.12400279 0.13013212 0.06024489
## CI 0.11537143 0.11643446 0.12400279 0.00000000 0.04780469 0.15838941
## TG 0.10498004 0.10810704 0.13013212 0.04780469 0.00000000 0.15012798
## TM 0.06321247 0.04798806 0.06024489 0.15838941 0.15012798 0.00000000
## 
## Pairwise 95% upper confidence limits between means
##            AG         AM         CF         CI         TG         TM
## AG 0.00000000 0.07675692 0.08148254 0.08379440 0.07699149 0.08157740
## AM 0.07675692 0.00000000 0.08216896 0.08273533 0.08113110 0.08009783
## CF 0.08148254 0.08216896 0.00000000 0.07738504 0.07915375 0.08217325
## CI 0.08379440 0.08273533 0.07738504 0.00000000 0.08263148 0.08728829
## TG 0.07699149 0.08113110 0.07915375 0.08263148 0.00000000 0.07862142
## TM 0.08157740 0.08009783 0.08217325 0.08728829 0.07862142 0.00000000
## 
## Pairwise effect sizes (Z) between means
##            AG         AM        CF        CI        TG        TM
## AG  0.0000000 -0.4761783 1.2593319 3.5878927 3.4659402 1.1375433
## AM -0.4761783  0.0000000 0.7376608 3.5504394 3.4251391 0.3995548
## CF  1.2593319  0.7376608 0.0000000 4.4056219 4.6313544 0.9170473
## CI  3.5878927  3.5504394 4.4056219 0.0000000 0.2730329 5.1622430
## TG  3.4659402  3.4251391 4.6313544 0.2730329 0.0000000 5.6782533
## TM  1.1375433  0.3995548 0.9170473 5.1622430 5.6782533 0.0000000
## 
## Pairwise P-values between means
##       AG    AM    CF    CI    TG    TM
## AG 1.000 0.594 0.121 0.005 0.006 0.130
## AM 0.594 1.000 0.202 0.007 0.005 0.262
## CF 0.121 0.202 1.000 0.003 0.001 0.183
## CI 0.005 0.007 0.003 1.000 0.325 0.001
## TG 0.006 0.005 0.001 0.325 1.000 0.001
## TM 0.130 0.262 0.183 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.011918019 0.012893856 0.015322630 0.013887004 0.008872023 0.011847448 
## 
## Pairwise distances between variances, plus statistics
##                  d   UCL (95%)           Z Pr > d
## AG:AM 9.758370e-04 0.007965799 -0.93063560  0.802
## AG:CF 3.404611e-03 0.007915377  0.06773634  0.404
## AG:CI 1.968985e-03 0.008222942 -0.55116824  0.641
## AG:TG 3.045996e-03 0.007965492 -0.10449792  0.474
## AG:TM 7.057081e-05 0.007725077 -1.29571098  0.988
## AM:CF 2.428774e-03 0.007865633 -0.38450554  0.576
## AM:CI 9.931479e-04 0.008102822 -0.91237916  0.811
## AM:TG 4.021833e-03 0.007847931  0.32864527  0.309
## AM:TM 1.046408e-03 0.008060068 -0.89107890  0.803
## CF:CI 1.435626e-03 0.008351726 -0.78896554  0.752
## CF:TG 6.450607e-03 0.008011795  1.23198800  0.121
## CF:TM 3.475182e-03 0.008158495  0.05377382  0.409
## CI:TG 5.014981e-03 0.008601174  0.58993706  0.258
## CI:TM 2.039556e-03 0.008189403 -0.49211358  0.613
## TG:TM 2.975425e-03 0.008579912 -0.15041950  0.469

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.011918019 0.012893856 0.015322630 0.013887004 0.008872023 0.011847448 
## 
## 
## Pairwise absolute differences between variances
##              AG           AM          CF           CI          TG
## AG 0.000000e+00 0.0009758370 0.003404611 0.0019689849 0.003045996
## AM 9.758370e-04 0.0000000000 0.002428774 0.0009931479 0.004021833
## CF 3.404611e-03 0.0024287738 0.000000000 0.0014356259 0.006450607
## CI 1.968985e-03 0.0009931479 0.001435626 0.0000000000 0.005014981
## TG 3.045996e-03 0.0040218328 0.006450607 0.0050149808 0.000000000
## TM 7.057081e-05 0.0010464078 0.003475182 0.0020395557 0.002975425
##              TM
## AG 7.057081e-05
## AM 1.046408e-03
## CF 3.475182e-03
## CI 2.039556e-03
## TG 2.975425e-03
## TM 0.000000e+00
## 
## 
## P-Values
##       AG    AM    CF    CI    TG    TM
## AG 1.000 0.802 0.404 0.641 0.474 0.988
## AM 0.802 1.000 0.576 0.811 0.309 0.803
## CF 0.404 0.576 1.000 0.752 0.121 0.409
## CI 0.641 0.811 0.752 1.000 0.258 0.613
## TG 0.474 0.309 0.121 0.258 1.000 0.469
## TM 0.988 0.803 0.409 0.613 0.469 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.2127472 0.1918373 0.6037024 0.2280649 0.2241217 0.1660466 
## 
## Pairwise absolute differences in path distances, plus statistics
##                 d UCL (95%)          Z Pr > d
## AG:AM 0.020909903 0.1924273 -0.8794736  0.797
## AG:CF 0.390955233 0.2126324  4.6960082  0.001
## AG:CI 0.015317794 0.2147393 -0.9611499  0.872
## AG:TG 0.011374573 0.2082791 -1.0774703  0.896
## AG:TM 0.046700541 0.1954840 -0.4957784  0.614
## AM:CF 0.411865136 0.2083110  5.1575623  0.001
## AM:CI 0.036227697 0.2058288 -0.6671303  0.698
## AM:TG 0.032284475 0.2053963 -0.7683800  0.753
## AM:TM 0.025790638 0.2105270 -0.8044252  0.771
## CF:CI 0.375637439 0.2088163  4.3203045  0.002
## CF:TG 0.379580660 0.2201008  4.3607967  0.001
## CF:TM 0.437655774 0.2158538  5.2246735  0.001
## CI:TG 0.003943221 0.2069388 -1.2063679  0.973
## CI:TM 0.062018335 0.2198893 -0.3513108  0.537
## TG:TM 0.058075114 0.2153191 -0.3266682  0.521

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.