Using RRPP
APPENDIX S2
This is a more detailed version of the section, “Analytical Demonstrations,” within the Methods in Ecology and Evolution article, “RRPP: An R package for fitting linear models to high-dimensional data using residual randomization,” by Michael L. Collyer and Dean C. Adams. The purpose of this vignette is to offer an opportunity for people reading that article to see code and results together. Some changes to code have been made, especially to enhance features that are difficult to introduce in the article, but the overall presentation should be basically the same. As the package has been updated, this document has also been slightly edited to make sure any functional changes have been illustrated.
(Note parallel processing is available on Unix systems for
lm.rrpp; see lm.rrpp help page for description
on its use. This is helpful for large data sets.)
Three data sets are provided in RRPP (with examples also
in help pages), and are used here as examples for various analyses:
Pupfish (Collyer et al., 2015)
PupfishHeads (Gilbert, 2016)
PlethMorph (Adams & Collyer, 2018)The first two datasets contain landmark-based geometric morphometric
data collected from museum samples of Pecos pupfish (Cyprinodon
pecosensis), representing body shape and cranial morphology,
respectively. Within both data sets, the $coords objects
are matrices of Procrustes residuals obtained from generalized
Procrustes analysis (GPA) of configurations of anatomical landmarks. For
the purposes of this example, it is sufficient to recognize that
Procrustes residuals embody a highly multivariate dataset representing
shape (see the R package, geomorph). The third data set contains
averaged linear measurements of 37 species of Plethodontid salamanders,
plus a covariance matrix based on a Brownian model of evolution, given
the phylogenetic relationship among the species (Adams & Collyer,
2018, for details).
Example: Pupfish Cranial Morphology and Mixed-Model ANOVA
In the first example we use a univariate dependent variable (head size, measured as the centroid size of the cranial landmark configuration; Bookstein, 1991) with a mixed-model design. The following code highlights the analytical steps, with results categorized in Fig. 1:
library(RRPP)
data("PupfishHeads")
PupfishHeads$logHeadSize <- log(PupfishHeads$headSize)
fit <- lm.rrpp(logHeadSize ~ sex + locality/year,
SS.type = "I", data = PupfishHeads,
print.progress = FALSE,
turbo = FALSE, verbose = TRUE)## Warning:
## This is not an error! It is a friendly warning.
##
## Because variables in the linear model are redundant,
## the linear model design has been truncated (via QR decomposition).
## Original X columns: 13
## Final X columns (rank): 9
## Check coefficients or degrees of freedom in ANOVA to see changes.
##
## Use suppressWarnings() to turn off these warnings.##
## Linear Model fit with lm.rrpp
##
## Number of observations: 268
## Number of dependent variables: 1
## Data space dimensions: 1
## Sums of Squares and Cross-products: Type I
## Number of permutations: 1000
##
## Full Model Analysis of Variance
##
## Df Residual Df SS Residual SS Rsq F
## sex + locality/year 8 259 2.718175 10.94887 0.1988854 8.037445
## Z (from F) Pr(>F)
## sex + locality/year 5.097397 0.001
##
##
## Redundancy Analysis (PCA on fitted values and residuals)
##
## Trace Proportion Rank
## Fitted 0.01018043 0.1988854 1
## Residuals 0.04100699 0.8011146 1
## Total 0.05118742 1.0000000 1
##
## Eigenvalues
##
## PC1
## Fitted 0.01018043
## Residuals 0.04100699
## Total 0.05118742##
## 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)
## sex 1 0.6091 0.60911 0.04457 14.4087 2.7924 0.002 **
## locality 1 0.4712 0.47121 0.03448 11.1466 2.6477 0.002 **
## locality:year 6 1.6379 0.27298 0.11984 6.4574 4.2051 0.001 ***
## Residuals 259 10.9489 0.04227 0.80111
## Total 267 13.6670
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Call: lm.rrpp(f1 = logHeadSize ~ sex + locality/year, turbo = FALSE,
## SS.type = "I", data = PupfishHeads, print.progress = FALSE,
## verbose = TRUE)Of important note, we choose to log-transform head size and include
it as a separate variable in the RRPP data frame,
PupfishHeads. We accomplished this via the code above –
rather than using log(headSize) - because downstream
functions like predict.lm.rrpp work better without
functions in the formula. ANOVA was performed using random distributions
of F-statistics to calculate z-scores and
P-values (but one could use alternative distributions - see the
RRPP help page). The S3 Generic functions
(summary, anova) return summaries that remind
the user how random data were generated, the type of SS, and
how z-scores were calculated. This particular ANOVA summary is
a default that fails to consider the year fish were sampled as a random
effect. A mixed-model ANOVA update can be performed by changing the
expected mean-square (MS) error estimates in each F
calculation:
##
## 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)
## sex 1 0.6091 0.60911 0.04457 14.4087 2.7924 0.002 **
## locality 1 0.4712 0.47121 0.03448 1.7262 0.7158 0.234
## locality:year 6 1.6379 0.27298 0.11984 6.4574 4.2051 0.001 ***
## Residuals 259 10.9489 0.04227 0.80111
## Total 267 13.6670
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Call: lm.rrpp(f1 = logHeadSize ~ sex + locality/year, turbo = FALSE,
## SS.type = "I", data = PupfishHeads, print.progress = FALSE,
## verbose = TRUE)This adjustment illustrates that the head size variation does not significantly differ between localities, with respect to the variation among sampling events. The anova function can also be used for multi-model inference.
fit.sex <- lm.rrpp(logHeadSize ~ sex,
data = PupfishHeads,
print.progress = FALSE)
fit.sex.loc<- lm.rrpp(logHeadSize ~ sex + locality,
data = PupfishHeads,
print.progress = FALSE)
fit.sex.loc.year<- lm.rrpp(logHeadSize ~ sex + locality/year,
data = PupfishHeads,
print.progress = FALSE)## Warning:
## This is not an error! It is a friendly warning.
##
## Because variables in the linear model are redundant,
## the linear model design has been truncated (via QR decomposition).
## Original X columns: 13
## Final X columns (rank): 9
## Check coefficients or degrees of freedom in ANOVA to see changes.
##
## Use suppressWarnings() to turn off these warnings.##
## 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
## logHeadSize ~ sex (Null) 266 1 13.058 0.000000
## logHeadSize ~ sex + locality 265 1 12.587 0.47121 0.47121 0.034478
## logHeadSize ~ sex + locality/year 259 7 10.949 2.10907 0.30130 0.154318
## Total 267 13.667
## F Z P Pr(>F)
## logHeadSize ~ sex (Null)
## logHeadSize ~ sex + locality 9.9207 2.5237 0.004
## logHeadSize ~ sex + locality/year 7.1273 4.5980 0.001
## TotalOne might wish to also look at individual model coefficients, and ascertain which have the largest effect:
##
## Linear Model fit with lm.rrpp
##
## Number of observations: 268
## Number of dependent variables: 1
## Data space dimensions: 1
## Sums of Squares and Cross-products: Type I
## Number of permutations: 1000
##
## Statistics (distances) of coefficients with 95 percent confidence intervals,
## effect sizes, and probabilities of exceeding observed values
## based on
## 1000 random permutations using RRPP
##
## d.obs UCL (95%) Zd Pr(>d)
## sexMale 0.09579938 0.05535808 2.6097241 0.003
## localitySH 0.09164189 0.05788327 2.5097854 0.004
## localitySH:year1998 0.06402169 0.13652635 0.4170950 0.355
## localityLake:year1999 0.35681331 0.13286852 3.7362990 0.001
## localitySH:year1999 0.01328613 0.13259002 -1.1032264 0.856
## localitySH:year2000 0.07165589 0.10621936 0.9403694 0.196
## localityLake:year2001 0.23035259 0.13911422 2.5419582 0.003
## localityLake:year2002 0.31050534 0.13479886 3.3397163 0.001This function produces a table much like summary.lm
output, but with bootstrap-generated confidence intervals of
coefficients.
It might be of interest to visualize model predictions for certain effects, holding constant other effects. For example, if we want to look at confidence intervals to compare male and female head sizes, holding constant the effects of locality and sampling period, we could do the following:
sizeDF <- data.frame(sex = c("Female", "Male"))
rownames(sizeDF) <- c("Female", "Male")
sizePreds <- predict(fit, sizeDF)## Warning:
## This is not an error! It is a friendly warning.
##
## Not all variables in model accounted for in newdata.
## Missing variables will be averaged from observed data for prediction.
##
## Use suppressWarnings to turn off these warnings.
The plots are perfectly amenable (e.g., point type and color, line
thickness, alternative labels, and additional text can be added or
adjusted with typical par arguments).
Finally, the SS type can be also toggled easily by refitting the model:
fit2 <- lm.rrpp(logHeadSize ~ sex + locality/year,
SS.type = "II", data = PupfishHeads, print.progress = FALSE)## Warning:
## This is not an error! It is a friendly warning.
##
## Because variables in the linear model are redundant,
## the linear model design has been truncated (via QR decomposition).
## Original X columns: 13
## Final X columns (rank): 9
## Check coefficients or degrees of freedom in ANOVA to see changes.
##
## Use suppressWarnings() to turn off these warnings.fit3 <- lm.rrpp(logHeadSize ~ sex + locality/year,
SS.type = "III", data = PupfishHeads, print.progress = FALSE)## Warning:
## This is not an error! It is a friendly warning.
##
## Because variables in the linear model are redundant,
## the linear model design has been truncated (via QR decomposition).
## Original X columns: 13
## Final X columns (rank): 9
## Check coefficients or degrees of freedom in ANOVA to see changes.
##
## Use suppressWarnings() to turn off these warnings.##
## 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)
## sex 1 0.6091 0.60911 0.04457 14.4087 2.7924 0.002 **
## locality 1 0.4712 0.47121 0.03448 11.1466 2.6477 0.002 **
## locality:year 6 1.6379 0.27298 0.11984 6.4574 4.2051 0.001 ***
## Residuals 259 10.9489 0.04227 0.80111
## Total 267 13.6670
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Call: lm.rrpp(f1 = logHeadSize ~ sex + locality/year, turbo = FALSE,
## SS.type = "I", data = PupfishHeads, print.progress = FALSE,
## verbose = TRUE)##
## 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 II
## Effect sizes (Z) based on F distributions
##
## Df SS MS Rsq F Z Pr(>F)
## sex 1 0.5554 0.55537 0.04064 13.1376 2.6737 0.002 **
## locality:year 6 1.6379 0.27298 0.11984 6.4574 4.2051 0.001 ***
## Residuals 259 10.9489 0.04227 0.80111
## Total 267 13.6670
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Call: lm.rrpp(f1 = logHeadSize ~ sex + locality/year, SS.type = "II",
## data = PupfishHeads, 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
## Sums of Squares and Cross-products: Type III
## Effect sizes (Z) based on F distributions
##
## Df SS MS Rsq F Z Pr(>F)
## sex 1 0.5554 0.55537 0.04064 13.1376 2.6737 0.002 **
## locality 1 0.0573 0.05733 0.00419 1.3562 0.7005 0.255
## locality:year 6 1.6379 0.27298 0.11984 6.4574 4.2051 0.001 ***
## Residuals 259 10.9489 0.04227 0.80111
## Total 267 13.6670
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Call: lm.rrpp(f1 = logHeadSize ~ sex + locality/year, SS.type = "III",
## data = PupfishHeads, print.progress = FALSE)Example: Pupfish Body Shape and High-Dimensional Data
In the second example, we highlight the RRPP ability to
efficiently handle large data computations. For this demonstration, a
54(n) × 112 (p) matrix of Procrustes residuals are the
data. In every one of the 1,000 random permutations, RRPP shuffles
residual vectors the same way for four different null models, estimates
coefficients for four different full models, estimates the SS
as the difference between residual SS (RSS) for four
null-full model comparisons, and calculates the total SS,
before calculating MS, R2, F,
Cohen’s f2, and Euclidean distances of coefficient
vectors across all 1,000 permutations. This process, plus packaging of
results, took approximately 0.5 seconds on a notebook computer, without
any parallel processing. In the second example, the initial steps are
quite the same as the first example:
data(Pupfish)
Pupfish$logSize <- log(Pupfish$CS)
fit <- lm.rrpp(coords ~ logSize + Sex*Pop, SS.type = "I",
data = Pupfish, print.progress = FALSE,
turbo = FALSE, verbose = TRUE)
summary(fit, formula = FALSE)##
## Linear Model fit with lm.rrpp
##
## Number of observations: 54
## Number of dependent variables: 112
## Data space dimensions: 53
## Sums of Squares and Cross-products: Type I
## Number of permutations: 1000
##
## Full Model Analysis of Variance
##
## Df Residual Df SS Residual SS Rsq F Z (from F) Pr(>F)
## fit 4 49 0.03224913 0.02408373 0.5724746 16.40327 7.783173 0.001
##
##
## Redundancy Analysis (PCA on fitted values and residuals)
##
## Trace Proportion Rank
## Fitted 0.0006084742 0.5724746 4
## Residuals 0.0004544100 0.4275254 49
## Total 0.0010628843 1.0000000 53
##
## Eigenvalues
##
## PC1 PC2 PC3 PC4 PC5
## Fitted 0.0003849870 0.0001764119 0.0000308588 0.0000162165
## Residuals 0.0001619045 0.0000683959 0.0000504732 0.0000375901 0.0000204706
## Total 0.0004644118 0.0002541398 0.0000722379 0.0000549026 0.0000491286
## PC6 PC7 PC8 PC9 PC10
## Fitted
## Residuals 0.0000158227 0.0000145861 0.0000120066 0.0000083853 0.0000074828
## Total 0.0000333998 0.0000227294 0.0000187381 0.0000127668 0.0000090893
## PC11 PC12 PC13 PC14 PC15
## Fitted
## Residuals 0.0000068974 0.0000062331 0.0000054188 0.0000045485 0.0000043871
## Total 0.0000088769 0.0000083661 0.0000064036 0.0000058026 0.0000051742
## PC16 PC17 PC18 PC19 PC20
## Fitted
## Residuals 0.0000032928 0.0000031713 0.0000027759 0.0000024382 0.0000023968
## Total 0.0000046349 0.0000039689 0.0000033385 0.0000029455 0.0000027293
## PC21 PC22 PC23 PC24 PC25
## Fitted
## Residuals 0.0000020800 0.0000016410 0.0000013892 0.0000012122 0.0000010647
## Total 0.0000024512 0.0000020110 0.0000019767 0.0000013932 0.0000012220
## PC26 PC27 PC28 PC29 PC30
## Fitted
## Residuals 0.0000009268 0.0000008974 0.0000008083 0.0000007337 0.0000006036
## Total 0.0000010286 0.0000010076 0.0000009490 0.0000008017 0.0000006669
## PC31 PC32 PC33 PC34 PC35
## Fitted
## Residuals 0.0000005856 0.0000004746 0.0000004236 0.0000004042 0.0000003381
## Total 0.0000006208 0.0000005973 0.0000005198 0.0000004633 0.0000004385
## PC36 PC37 PC38 PC39 PC40
## Fitted
## Residuals 0.0000003127 0.0000002573 0.0000002486 0.0000002226 0.0000002059
## Total 0.0000003599 0.0000003568 0.0000003050 0.0000002580 0.0000002252
## PC41 PC42 PC43 PC44 PC45
## Fitted
## Residuals 0.0000001783 0.0000001509 0.0000001282 0.0000001152 0.0000000896
## Total 0.0000002162 0.0000001941 0.0000001842 0.0000001527 0.0000001461
## PC46 PC47 PC48 PC49 PC50
## Fitted
## Residuals 0.0000000737 0.0000000708 0.0000000488 0.0000000468
## Total 0.0000001100 0.0000000963 0.0000000850 0.0000000697 0.0000000563
## PC51 PC52 PC53
## Fitted
## Residuals
## Total 0.0000000525 0.0000000447 0.0000000395##
## 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)
## logSize 1 0.014019 0.0140193 0.24886 28.5232 4.9982 0.001 ***
## Sex 1 0.007615 0.0076151 0.13518 15.4933 4.6798 0.001 ***
## Pop 1 0.007661 0.0076606 0.13599 15.5860 4.3464 0.001 ***
## Sex:Pop 1 0.002954 0.0029542 0.05244 6.0105 3.4054 0.002 **
## Residuals 49 0.024084 0.0004915 0.42753
## Total 53 0.056333
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Call: lm.rrpp(f1 = coords ~ logSize + Sex * Pop, turbo = FALSE, SS.type = "I",
## data = Pupfish, print.progress = FALSE, verbose = TRUE)##
## Linear Model fit with lm.rrpp
##
## Number of observations: 54
## Number of dependent variables: 112
## Data space dimensions: 53
## Sums of Squares and Cross-products: Type I
## Number of permutations: 1000
##
## Statistics (distances) of coefficients with 95 percent confidence intervals,
## effect sizes, and probabilities of exceeding observed values
## based on
## 1000 random permutations using RRPP
##
## d.obs UCL (95%) Zd Pr(>d)
## logSize 0.11101234 0.04871922 4.148348 0.001
## SexM 0.02755503 0.01300605 4.120183 0.001
## PopSinkhole 0.02895656 0.01256411 4.121852 0.001
## SexM:PopSinkhole 0.03002431 0.01772916 3.469951 0.002ANOVA results reveal that after accounting for body size allometry, not only are there significant inter-population differences in body shape and sexual dimorphism in body shape, but sexual dimorphism also significantly varies between the two populations. A fuller evaluation of these results is available in Collyer et al. (2015). It is worth taking a moment to realize why this analysis is a valuable alternative to a parametric M-ANOVA. The following code attempts to perform a parametric M-ANOVA:
fit$LM$data$coords <- Pupfish$coords
fit.par <- lm(fit$call$f1, data = fit$LM$data)
all.equal(fit$LM$coefficients, fit.par$coefficients)## [1] TRUE## Error in summary.manova(manova(fit.par)): residuals have rank 49 < 112Although both functions return the same coefficients, the error
summarizing the attempted parametric M-ANOVA clearly indicates the
limitation of having residual degrees of freedom (rank) lower than the
number of variables. Returning to the lm.rrpp fit, which
does not suffer this problem, we can look at the precision of group mean
estimation, accounting for allometric shape variation, by doing the
following:
shapeDF <- expand.grid(Sex = levels(Pupfish$Sex), Pop = levels(Pupfish$Pop))
rownames(shapeDF) <- paste(shapeDF$Sex, shapeDF$Pop, sep = ".")
shapePreds <- predict(fit, shapeDF, confidence = 0.95)## Warning:
## This is not an error! It is a friendly warning.
##
## Not all variables in model accounted for in newdata.
## Missing variables will be averaged from observed data for prediction.
##
## Use suppressWarnings to turn off these warnings.
plot(shapePreds, PC = TRUE, ellipse = TRUE,
pch = 19, col = 1:NROW(shapeDF)) # with added par arguments
These plots differ as there is a rotational difference between the
covariance matrices estimated with 4 predicted and 54 fitted values.
Additionally, the former illustrates prediction precision and the latter
sample dispersion. Both functions allow passing par
arguments to the plot as well as saving plot data for more advanced
plotting. The following is a regression-type plot:
plot(fit, type = "regression", reg.type = "PredLine",
predictor = Pupfish$logSize, pch=19,
col = as.numeric(groups))
The function, pairwise, can be used to test pairwise differences between least-squares means with:
PWT <- pairwise(fit, groups = interaction(Pupfish$Sex, Pupfish$Pop))
summary(PWT, confidence = 0.95)##
## Pairwise comparisons
##
## Groups: F.Marsh M.Marsh F.Sinkhole M.Sinkhole
##
## 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
## F.Marsh:M.Marsh 0.03579214 0.03278645 2.2202630 0.015
## F.Marsh:F.Sinkhole 0.03639758 0.03601481 1.7807393 0.044
## F.Marsh:M.Sinkhole 0.03809715 0.04432235 -0.6390591 0.735
## M.Marsh:F.Sinkhole 0.03681731 0.04779814 -0.2465676 0.601
## M.Marsh:M.Sinkhole 0.02694896 0.03646576 -0.9484513 0.830
## F.Sinkhole:M.Sinkhole 0.01939739 0.03295865 -1.7845209 0.962Much like the tukeyHSD function in the R stats package,
pairwise will generate tables with confidence intervals and
P-values for the pairwise statistic, Euclidean distance between
least-squares means. This function could also be used for pairwise
comparison of slopes in analysis of covariance (ANCOVA) designs.
fit2 <- lm.rrpp(coords ~ logSize * Sex * Pop, SS.type = "I",
data = Pupfish, print.progress = FALSE, iter = 999)
summary(fit2, formula = FALSE)##
## Linear Model fit with lm.rrpp
##
## Number of observations: 54
## Number of dependent variables: 112
## Data space dimensions: 53
## Sums of Squares and Cross-products: Type I
## Number of permutations: 1000
##
## Full Model Analysis of Variance
##
## Df Residual Df SS Residual SS Rsq F Z (from F) Pr(>F)
## fit2 7 46 0.03425362 0.02207925 0.6080574 10.19488 9.289796 0.001
##
##
## Redundancy Analysis (PCA on fitted values and residuals)
##
## Trace Proportion Rank
## Fitted 0.0006462947 0.6080575 7
## Residuals 0.0004165896 0.3919426 46
## Total 0.0010628843 1.0000000 53
##
## Eigenvalues
##
## PC1 PC2 PC3 PC4 PC5
## Fitted 0.0003957338 0.0001850501 0.0000318297 0.0000187072 0.0000081143
## Residuals 0.0001395372 0.0000679693 0.0000501429 0.0000336509 0.0000189840
## Total 0.0004644118 0.0002541398 0.0000722379 0.0000549026 0.0000491286
## PC6 PC7 PC8 PC9 PC10
## Fitted 0.0000038627 0.0000029970
## Residuals 0.0000153088 0.0000133860 0.0000109149 0.0000081753 0.0000072355
## Total 0.0000333998 0.0000227294 0.0000187381 0.0000127668 0.0000090893
## PC11 PC12 PC13 PC14 PC15
## Fitted
## Residuals 0.0000066482 0.0000053575 0.0000051028 0.0000040219 0.0000034070
## Total 0.0000088769 0.0000083661 0.0000064036 0.0000058026 0.0000051742
## PC16 PC17 PC18 PC19 PC20
## Fitted
## Residuals 0.0000030820 0.0000029693 0.0000025181 0.0000024262 0.0000019439
## Total 0.0000046349 0.0000039689 0.0000033385 0.0000029455 0.0000027293
## PC21 PC22 PC23 PC24 PC25
## Fitted
## Residuals 0.0000018143 0.0000016024 0.0000013239 0.0000011070 0.0000009227
## Total 0.0000024512 0.0000020110 0.0000019767 0.0000013932 0.0000012220
## PC26 PC27 PC28 PC29 PC30
## Fitted
## Residuals 0.0000008477 0.0000008014 0.0000007148 0.0000006218 0.0000005587
## Total 0.0000010286 0.0000010076 0.0000009490 0.0000008017 0.0000006669
## PC31 PC32 PC33 PC34 PC35
## Fitted
## Residuals 0.0000005125 0.0000004234 0.0000003470 0.0000003184 0.0000002807
## Total 0.0000006208 0.0000005973 0.0000005198 0.0000004633 0.0000004385
## PC36 PC37 PC38 PC39 PC40
## Fitted
## Residuals 0.0000002721 0.0000002379 0.0000002230 0.0000001699 0.0000001561
## Total 0.0000003599 0.0000003568 0.0000003050 0.0000002580 0.0000002252
## PC41 PC42 PC43 PC44 PC45
## Fitted
## Residuals 0.0000001439 0.0000001208 0.0000000942 0.0000000763 0.0000000673
## Total 0.0000002162 0.0000001941 0.0000001842 0.0000001527 0.0000001461
## PC46 PC47 PC48 PC49 PC50
## Fitted
## Residuals 0.0000000496
## Total 0.0000001100 0.0000000963 0.0000000850 0.0000000697 0.0000000563
## PC51 PC52 PC53
## Fitted
## Residuals
## Total 0.0000000525 0.0000000447 0.0000000395##
## 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 ~ logSize + Sex * Pop (Null) 49 1 0.024084
## coords ~ logSize * Sex * Pop 46 3 0.022079 0.0020045 0.00066816
## Total 53 0.056333
## Rsq F Z P Pr(>F)
## coords ~ logSize + Sex * Pop (Null) 0.000000
## coords ~ logSize * Sex * Pop 0.035583 1.3921 1.1451 0.136
## TotalPW2 <- pairwise(fit2, fit.null = fit, groups = groups,
covariate = Pupfish$logSize, print.progress = FALSE)
PW2##
## Pairwise comparisons
##
## Groups: F.Marsh M.Marsh F.Sinkhole M.Sinkhole
##
## RRPP: 1000 permutations##
## Pairwise comparisons
##
## Groups: F.Marsh M.Marsh F.Sinkhole M.Sinkhole
##
## RRPP: 1000 permutations
##
## Slopes (vectors of variate change per one unit of covariate
## change, by group):
## Vectors hidden (use show.vectors = TRUE to view)
##
## Pairwise distances between slope vector
## (end-points), plus statistics
## d UCL (95%) Z Pr > d
## F.Marsh:M.Marsh 0.08279940 0.1640452 -0.9689342 0.829
## F.Marsh:F.Sinkhole 0.07848121 0.1621176 -0.9763904 0.838
## F.Marsh:M.Sinkhole 0.12322729 0.1610484 0.7139190 0.236
## M.Marsh:F.Sinkhole 0.06418279 0.1095190 -0.3981426 0.652
## M.Marsh:M.Sinkhole 0.10863856 0.1132783 1.5445168 0.063
## F.Sinkhole:M.Sinkhole 0.11706696 0.1001323 2.2066562 0.015##
## Pairwise comparisons
##
## Groups: F.Marsh M.Marsh F.Sinkhole M.Sinkhole
##
## RRPP: 1000 permutations
##
## Slopes (vectors of variate change per one unit of covariate
## change, by group):
## Vectors hidden (use show.vectors = TRUE to view)
##
## Pairwise distances between slope vector (end-points
## F.Marsh M.Marsh F.Sinkhole M.Sinkhole
## F.Marsh 0.00000000 0.08279940 0.07848121 0.1232273
## M.Marsh 0.08279940 0.00000000 0.06418279 0.1086386
## F.Sinkhole 0.07848121 0.06418279 0.00000000 0.1170670
## M.Sinkhole 0.12322729 0.10863856 0.11706696 0.0000000
##
## Pairwise 95% Upper confidence limits between slopes
## F.Marsh M.Marsh F.Sinkhole M.Sinkhole
## F.Marsh 0.0000000 0.1640452 0.1621176 0.1610484
## M.Marsh 0.1640452 0.0000000 0.1095190 0.1132783
## F.Sinkhole 0.1621176 0.1095190 0.0000000 0.1001323
## M.Sinkhole 0.1610484 0.1132783 0.1001323 0.0000000
##
## Pairwise effect sizes (Z) between slopes
## F.Marsh M.Marsh F.Sinkhole M.Sinkhole
## F.Marsh 0.0000000 -0.9689342 -0.9763904 0.713919
## M.Marsh -0.9689342 0.0000000 -0.3981426 1.544517
## F.Sinkhole -0.9763904 -0.3981426 0.0000000 2.206656
## M.Sinkhole 0.7139190 1.5445168 2.2066562 0.000000
##
## Pairwise P-values between slopes
## F.Marsh M.Marsh F.Sinkhole M.Sinkhole
## F.Marsh 1.000 0.829 0.838 0.236
## M.Marsh 0.829 1.000 0.652 0.063
## F.Sinkhole 0.838 0.652 1.000 0.015
## M.Sinkhole 0.236 0.063 0.015 1.000summary(PW2, confidence = 0.95,
test.type = "VC",
angle.type = "deg") # correlation between slope vectors (and angles)##
## Pairwise comparisons
##
## Groups: F.Marsh M.Marsh F.Sinkhole M.Sinkhole
##
## RRPP: 1000 permutations
##
## Slopes (vectors of variate change per one unit of covariate
## change, by group):
## Vectors hidden (use show.vectors = TRUE to view)
##
## Pairwise statistics based on slopes vector correlations (r)
## and angles, acos(r)
## The null hypothesis is that r = 1 (parallel vectors).
## This null hypothesis is better treated as the angle
## between vectors = 0
## r angle UCL (95%) Z Pr > angle
## F.Marsh:M.Marsh 0.6139591 52.12367 83.58445 -0.13196270 0.549
## F.Marsh:F.Sinkhole 0.6318267 50.81498 83.57546 -0.07310658 0.527
## F.Marsh:M.Sinkhole 0.4461018 63.50614 81.89039 0.66791084 0.273
## M.Marsh:F.Sinkhole 0.7175129 44.15048 63.60905 0.24275646 0.391
## M.Marsh:M.Sinkhole 0.5629975 55.73665 65.19951 1.04510764 0.154
## F.Sinkhole:M.Sinkhole 0.4627734 62.43378 60.22153 1.77232308 0.037Because the Procrustes residuals are projected into a Euclidean
tangent space (see geomorph function, gpagen; Adams et al., 2017), this
analysis could be performed with an object of class dist
(values from lower half of a distance matrix) representing the
inter-specimen shape (Euclidean) distances, using the following
code:
D <- dist(Pupfish$coords) # inter-observation Euclidean distances
Pupfish$D <- D
fitD <- lm.rrpp(D ~ logSize + Sex*Pop, SS.type = "I",
data = Pupfish, print.progress = FALSE)
summary(fitD)##
## Linear Model fit with lm.rrpp
##
## Number of observations: 54
## Number of dependent variables: 53
## Data space dimensions: 53
## Sums of Squares and Cross-products: Type I
## Number of permutations: 1000
##
## Full Model Analysis of Variance
##
## Df Residual Df SS Residual SS Rsq F
## logSize + Sex * Pop 4 49 0.03224913 0.02408373 0.5724746 16.40327
## Z (from F) Pr(>F)
## logSize + Sex * Pop 7.783173 0.001
##
##
## Redundancy Analysis (PCA on fitted values and residuals)
##
## Trace Proportion Rank
## Fitted 0.0006084742 0.5724746 4
## Residuals 0.0004544100 0.4275254 49
## Total 0.0010628843 1.0000000 53
##
## Eigenvalues
##
## PC1 PC2 PC3 PC4 PC5
## Fitted 0.0003849870 0.0001764119 0.0000308588 0.0000162165
## Residuals 0.0001619045 0.0000683959 0.0000504732 0.0000375901 0.0000204706
## Total 0.0004644118 0.0002541398 0.0000722379 0.0000549026 0.0000491286
## PC6 PC7 PC8 PC9 PC10
## Fitted
## Residuals 0.0000158227 0.0000145861 0.0000120066 0.0000083853 0.0000074828
## Total 0.0000333998 0.0000227294 0.0000187381 0.0000127668 0.0000090893
## PC11 PC12 PC13 PC14 PC15
## Fitted
## Residuals 0.0000068974 0.0000062331 0.0000054188 0.0000045485 0.0000043871
## Total 0.0000088769 0.0000083661 0.0000064036 0.0000058026 0.0000051742
## PC16 PC17 PC18 PC19 PC20
## Fitted
## Residuals 0.0000032928 0.0000031713 0.0000027759 0.0000024382 0.0000023968
## Total 0.0000046349 0.0000039689 0.0000033385 0.0000029455 0.0000027293
## PC21 PC22 PC23 PC24 PC25
## Fitted
## Residuals 0.0000020800 0.0000016410 0.0000013892 0.0000012122 0.0000010647
## Total 0.0000024512 0.0000020110 0.0000019767 0.0000013932 0.0000012220
## PC26 PC27 PC28 PC29 PC30
## Fitted
## Residuals 0.0000009268 0.0000008974 0.0000008083 0.0000007337 0.0000006036
## Total 0.0000010286 0.0000010076 0.0000009490 0.0000008017 0.0000006669
## PC31 PC32 PC33 PC34 PC35
## Fitted
## Residuals 0.0000005856 0.0000004746 0.0000004236 0.0000004042 0.0000003381
## Total 0.0000006208 0.0000005973 0.0000005198 0.0000004633 0.0000004385
## PC36 PC37 PC38 PC39 PC40
## Fitted
## Residuals 0.0000003127 0.0000002573 0.0000002486 0.0000002226 0.0000002059
## Total 0.0000003599 0.0000003568 0.0000003050 0.0000002580 0.0000002252
## PC41 PC42 PC43 PC44 PC45
## Fitted
## Residuals 0.0000001783 0.0000001509 0.0000001282 0.0000001152 0.0000000896
## Total 0.0000002162 0.0000001941 0.0000001842 0.0000001527 0.0000001461
## PC46 PC47 PC48 PC49 PC50
## Fitted
## Residuals 0.0000000737 0.0000000708 0.0000000488 0.0000000468
## Total 0.0000001100 0.0000000963 0.0000000850 0.0000000697 0.0000000563
## PC51 PC52 PC53
## Fitted
## Residuals
## Total 0.0000000525 0.0000000447 0.0000000395##
## Linear Model fit with lm.rrpp
##
## Number of observations: 54
## Number of dependent variables: 112
## Data space dimensions: 53
## Sums of Squares and Cross-products: Type I
## Number of permutations: 1000
##
## Full Model Analysis of Variance
##
## Df Residual Df SS Residual SS Rsq F
## logSize + Sex * Pop 4 49 0.03224913 0.02408373 0.5724746 16.40327
## Z (from F) Pr(>F)
## logSize + Sex * Pop 7.783173 0.001
##
##
## Redundancy Analysis (PCA on fitted values and residuals)
##
## Trace Proportion Rank
## Fitted 0.0006084742 0.5724746 4
## Residuals 0.0004544100 0.4275254 49
## Total 0.0010628843 1.0000000 53
##
## Eigenvalues
##
## PC1 PC2 PC3 PC4 PC5
## Fitted 0.0003849870 0.0001764119 0.0000308588 0.0000162165
## Residuals 0.0001619045 0.0000683959 0.0000504732 0.0000375901 0.0000204706
## Total 0.0004644118 0.0002541398 0.0000722379 0.0000549026 0.0000491286
## PC6 PC7 PC8 PC9 PC10
## Fitted
## Residuals 0.0000158227 0.0000145861 0.0000120066 0.0000083853 0.0000074828
## Total 0.0000333998 0.0000227294 0.0000187381 0.0000127668 0.0000090893
## PC11 PC12 PC13 PC14 PC15
## Fitted
## Residuals 0.0000068974 0.0000062331 0.0000054188 0.0000045485 0.0000043871
## Total 0.0000088769 0.0000083661 0.0000064036 0.0000058026 0.0000051742
## PC16 PC17 PC18 PC19 PC20
## Fitted
## Residuals 0.0000032928 0.0000031713 0.0000027759 0.0000024382 0.0000023968
## Total 0.0000046349 0.0000039689 0.0000033385 0.0000029455 0.0000027293
## PC21 PC22 PC23 PC24 PC25
## Fitted
## Residuals 0.0000020800 0.0000016410 0.0000013892 0.0000012122 0.0000010647
## Total 0.0000024512 0.0000020110 0.0000019767 0.0000013932 0.0000012220
## PC26 PC27 PC28 PC29 PC30
## Fitted
## Residuals 0.0000009268 0.0000008974 0.0000008083 0.0000007337 0.0000006036
## Total 0.0000010286 0.0000010076 0.0000009490 0.0000008017 0.0000006669
## PC31 PC32 PC33 PC34 PC35
## Fitted
## Residuals 0.0000005856 0.0000004746 0.0000004236 0.0000004042 0.0000003381
## Total 0.0000006208 0.0000005973 0.0000005198 0.0000004633 0.0000004385
## PC36 PC37 PC38 PC39 PC40
## Fitted
## Residuals 0.0000003127 0.0000002573 0.0000002486 0.0000002226 0.0000002059
## Total 0.0000003599 0.0000003568 0.0000003050 0.0000002580 0.0000002252
## PC41 PC42 PC43 PC44 PC45
## Fitted
## Residuals 0.0000001783 0.0000001509 0.0000001282 0.0000001152 0.0000000896
## Total 0.0000002162 0.0000001941 0.0000001842 0.0000001527 0.0000001461
## PC46 PC47 PC48 PC49 PC50
## Fitted
## Residuals 0.0000000737 0.0000000708 0.0000000488 0.0000000468
## Total 0.0000001100 0.0000000963 0.0000000850 0.0000000697 0.0000000563
## PC51 PC52 PC53
## Fitted
## Residuals
## Total 0.0000000525 0.0000000447 0.0000000395##
## 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)
## logSize 1 0.014019 0.0140193 0.24886 28.5232 4.9982 0.001 ***
## Sex 1 0.007615 0.0076151 0.13518 15.4933 4.6798 0.001 ***
## Pop 1 0.007661 0.0076606 0.13599 15.5860 4.3464 0.001 ***
## Sex:Pop 1 0.002954 0.0029542 0.05244 6.0105 3.4054 0.002 **
## Residuals 49 0.024084 0.0004915 0.42753
## Total 53 0.056333
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Call: lm.rrpp(f1 = D ~ logSize + Sex * Pop, SS.type = "I", data = Pupfish,
## 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
## Sums of Squares and Cross-products: Type I
## Effect sizes (Z) based on F distributions
##
## Df SS MS Rsq F Z Pr(>F)
## logSize 1 0.014019 0.0140193 0.24886 28.5232 4.9982 0.001 ***
## Sex 1 0.007615 0.0076151 0.13518 15.4933 4.6798 0.001 ***
## Pop 1 0.007661 0.0076606 0.13599 15.5860 4.3464 0.001 ***
## Sex:Pop 1 0.002954 0.0029542 0.05244 6.0105 3.4054 0.002 **
## Residuals 49 0.024084 0.0004915 0.42753
## Total 53 0.056333
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Call: lm.rrpp(f1 = coords ~ logSize + Sex * Pop, turbo = FALSE, SS.type = "I",
## data = Pupfish, print.progress = FALSE, verbose = TRUE)The ANOVA results with either method are exactly the same.
Example: Plethodontid Morphology, Phylogenetics, and GLS Estimation
In the third example, we highlight GLS estimation. The following code
creates two lm.rrpp fits using OLS and GLS, respectively,
and evaluates them as in previous examples:
data(PlethMorph)
fitOLS <- lm.rrpp(TailLength ~ SVL,
data = PlethMorph,
print.progress = FALSE,
turbo = FALSE, verbose = TRUE)
fitGLS <- lm.rrpp(TailLength ~ SVL,
data = PlethMorph,
Cov = PlethMorph$PhyCov,
print.progress = FALSE,
turbo = FALSE, verbose = TRUE)
anova(fitOLS)##
## 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)
## SVL 1 3707.9 3707.9 0.76849 116.18 5.7522 0.001 ***
## Residuals 35 1117.0 31.9 0.23151
## Total 36 4824.9
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Call: lm.rrpp(f1 = TailLength ~ SVL, turbo = FALSE, data = PlethMorph,
## print.progress = FALSE, verbose = TRUE)##
## Analysis of Variance, using Residual Randomization
## Permutation procedure: Randomization of null model residuals
## Number of permutations: 1000
## Estimation method: Generalized Least-Squares (via OLS projection)
## Sums of Squares and Cross-products: Type I
## Effect sizes (Z) based on F distributions
##
## Df SS MS Rsq F Z Pr(>F)
## SVL 1 189.34 189.343 0.34259 18.239 3.2047 0.001 ***
## Residuals 35 363.33 10.381 0.65741
## Total 36 552.68
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Call: lm.rrpp(f1 = TailLength ~ SVL, turbo = FALSE, data = PlethMorph,
## Cov = PlethMorph$PhyCov, print.progress = FALSE, verbose = TRUE)##
## Linear Model fit with lm.rrpp
##
## Number of observations: 37
## Number of dependent variables: 1
## Data space dimensions: 1
## Sums of Squares and Cross-products: Type I
## Number of permutations: 1000
##
## Statistics (distances) of coefficients with 95 percent confidence intervals,
## effect sizes, and probabilities of exceeding observed values
## based on
## 1000 random permutations using RRPP
##
## d.obs UCL (95%) Zd Pr(>d)
## SVL 0.9681511 0.3597148 3.907365 0.001##
## Linear Model fit with lm.rrpp
##
## Number of observations: 37
## Number of dependent variables: 1
## Data space dimensions: 1
## Sums of Squares and Cross-products: Type I
## Number of permutations: 1000
##
## Statistics (distances) of coefficients with 95 percent confidence intervals,
## effect sizes, and probabilities of exceeding observed values
## based on
## 1000 random permutations using RRPP
##
## d.obs UCL (95%) Zd Pr(>d)
## SVL 0.8307223 0.4503532 2.894442 0.001Although analyses on either model fit indicate a significant relationship between tail length and snout-to-vent length (SVL), the GLS coefficients test and ANOVA show how phylogenetic auto-correlation among species augments the OLS-estimated relationship. The following is a multivariate example:
Y <- as.matrix(cbind(PlethMorph$TailLength,
PlethMorph$HeadLength,
PlethMorph$TailLength,
PlethMorph$Snout.eye,
PlethMorph$BodyWidth,
PlethMorph$Forelimb,
PlethMorph$Hindlimb))
PlethMorph$Y <- Y
fitOLSm <- lm.rrpp(Y ~ SVL, data = PlethMorph,
print.progress = FALSE,
turbo = FALSE, verbose = TRUE)
fitGLSm <- lm.rrpp(Y ~ SVL, data = PlethMorph,
Cov = PlethMorph$PhyCov,
print.progress = FALSE,
turbo = FALSE, verbose = TRUE)
anova(fitOLSm)##
## 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)
## SVL 1 8135.3 8135.3 0.76867 116.3 5.0898 0.001 ***
## Residuals 35 2448.3 70.0 0.23133
## Total 36 10583.7
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Call: lm.rrpp(f1 = Y ~ SVL, turbo = FALSE, data = PlethMorph, print.progress = FALSE,
## verbose = TRUE)##
## Analysis of Variance, using Residual Randomization
## Permutation procedure: Randomization of null model residuals
## Number of permutations: 1000
## Estimation method: Generalized Least-Squares (via OLS projection)
## Sums of Squares and Cross-products: Type I
## Effect sizes (Z) based on F distributions
##
## Df SS MS Rsq F Z Pr(>F)
## SVL 1 395.35 395.35 0.34715 18.611 3.0419 0.001 ***
## Residuals 35 743.49 21.24 0.65285
## Total 36 1138.84
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Call: lm.rrpp(f1 = Y ~ SVL, turbo = FALSE, data = PlethMorph, Cov = PlethMorph$PhyCov,
## print.progress = FALSE, verbose = TRUE)sizeDF <- data.frame(SVL = sort(PlethMorph$SVL))
plot(predict(fitOLSm, sizeDF), PC= TRUE) # Correlated error
Analytical Summary
On the surface, these three examples and their analyses should seem
intuitive to any user of R who has used the lm function
plus its associated S3 generics
(coef, predict, resid, fitted, summary, and
anova), all of which can be used on lm.rrpp
model fits. The functions, pairwise (not an S3 generic) and
anova, also allow pairwise comparisons of least-squares
means or slopes and multi-model inferences, respectively. Advanced users
will recognize, however, much more extensive usable results for adaptive
programming. The output from a lm.rrpp fit is arranged
hierarchically, i.e.:
## $names
## [1] "call" "LM" "ANOVA" "PermInfo" "turbo" "Models" "verbose"
## [8] "subTest"
##
## $class
## [1] "lm.rrpp"Within the $LM partition, all attributes of the lm
function are found, in addition to coefficients for every random
permutation. Within the $ANOVA partition, the SS type, plus
SS, MS, R2,
F, and Cohen’s f2 for all
permutations, as well as effect sizes estimated for each of these are
provided. Within the $PermInfo partition, the number of
permutations, type (RRPP or randomization of “full” data values, FRPP),
and sampling frame in every permutation (schedule) are provided. Thus,
lm.rrpp is the workhorse that makes all downstream analysis
efficient.
References
Adams, D. C., & Collyer, M. L. (2018). Multivariate phylogenetic comparative methods: Evaluations, comparisons, and recommendations. Systematic Biology, 67, 14–31.
Adams, D. C., Collyer, M. L., Kaliontzopoulou, A., & Sherratt, E. (2017). Geomorph: Software for geometric morphometric analyses. R package version 3.0.6.
Adams, D. C., Rohlf, F. J., & Slice, D. E. 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: Cambridge University Press.
Collyer, M. L., Sekora, D. J., & Adams, D. C. (2015). A method for analysis of phenotypic change for phenotypes described by high-dimensional data. Heredity, 115, 357–365.
Gilbert, M. C. (2016). Impacts of habitat fragmentation on the cranial morphology of a threatened desert fish (Cyprinodon pecosensis). Masters thesis. Western Kentucky University, Bowling Green, KY, USA.