Characterizing patterns of integration and modularity are hugely important in morphometrics. At their core, these methods seek to summarize patterns of trait covariation, as found in the trait covariance matrix: \(\hat{\mathbf\Sigma} = \mathbf{Z^T}\mathbf{Z}/ (n-1)\). There are three general classes of methods:
`Geomorph implements approaches addressing all three of these (See lecture for detailed description of methods).
library(geomorph)
## Loading required package: RRPP
## Loading required package: rgl
## Loading required package: Matrix
Here we use an effect size based on \(V_{rel}\) to quantify and compare overall levels of integration in a set of traits.
data("plethodon")
Y.gpa <- gpagen(plethodon$land, print.progress = FALSE)
#Separate data by species
coords.gp <- coords.subset(Y.gpa$coords, plethodon$species)
#Z_Vrel by species
Vrel.gp <- Map(function(x) integration.Vrel(x), coords.gp)
compare.ZVrel(Vrel.gp$Jord, Vrel.gp$Teyah)
##
## Effect sizes
##
## Vrel.gp$Jord Vrel.gp$Teyah
## -0.2978931 -0.2642648
##
## Effect sizes for pairwise differences in rel.eig effect size
##
## Vrel.gp$Jord Vrel.gp$Teyah
## Vrel.gp$Jord 0.00000000 0.09804249
## Vrel.gp$Teyah 0.09804249 0.00000000
##
## P-values
##
## Vrel.gp$Jord Vrel.gp$Teyah
## Vrel.gp$Jord 1.0000000 0.9218986
## Vrel.gp$Teyah 0.9218986 1.0000000
A method for evaluating integration across spatial scales is also implemented
globalIntegration(Y.gpa$coords) #data are not spatially integrated
## BEval
## -0.6369198
Here, one may use PLS to evaluate the covariation between subsets of traits:
data(pupfish) # GPA previously performed
group <- factor(paste(pupfish$Pop, pupfish$Sex, sep = "."))
# Subset 3D array by group, returning a list of 3D arrays
tail.LM <- c(1:3, 5:9, 18:38)
head.LM <- (1:56)[-tail.LM]
tail.coords <- pupfish$coords[tail.LM,,]
head.coords <- pupfish$coords[head.LM,,]
IT <- integration.test(tail.coords, head.coords, print.progress = F)
##
## No names for A2. Data are assumed to be ordered as in A.
summary(IT)
##
## Call:
## integration.test(A = tail.coords, A2 = head.coords, print.progress = F)
##
##
##
## r-PLS: 0.917
##
## Effect Size (Z): 5.4983
##
## P-value: 0.001
##
## Based on 1000 random permutations
plot(IT)
One could also use a partition of landmarks (NOTE identical \(r_{PLS}\):
land.gp <- rep(1,56); land.gp[tail.LM] <- 2
integration.test(pupfish$coords, partition.gp=land.gp, print.progress = FALSE)
##
## Call:
## integration.test(A = pupfish$coords, partition.gp = land.gp,
## print.progress = FALSE)
##
##
##
## r-PLS: 0.917
##
## Effect Size (Z): 5.4019
##
## P-value: 0.001
##
## Based on 1000 random permutations
And the statistical results are the same as when using two.b.pls
two.b.pls(tail.coords, head.coords, print.progress = FALSE)
## Data in either A1 or A2 do not have names. It is assumed data in both A1 and A2 are ordered the same.
##
## Call:
## two.b.pls(A1 = tail.coords, A2 = head.coords, print.progress = FALSE)
##
##
##
## r-PLS: 0.917
##
## Effect Size (Z): 5.4983
##
## P-value: 0.001
##
## Based on 1000 random permutations
picknplot.shape
example is provided in lab script.
One may also evaluate whether the strength of integration similar in two or more groups. As above this is accomplished using effect sizes, this time based on \(r_{PLS}\).
tail.coords.gp <- coords.subset(tail.coords, group)
head.coords.gp <- coords.subset(head.coords, group)
# Obtain Integration for groups
integ.tests <- Map(function(x,y) integration.test(x, y, iter=499,
print.progress = FALSE), head.coords.gp, tail.coords.gp)
compare.pls(integ.tests)
##
## Effect sizes
##
## Marsh.F Marsh.M Sinkhole.F Sinkhole.M
## 3.1410345 1.5130737 0.7429882 2.7156556
##
## Effect sizes for pairwise differences in PLS effect size
##
## Marsh.F Marsh.M Sinkhole.F Sinkhole.M
## Marsh.F 0.0000000 0.6586628 1.9027247 0.7518767
## Marsh.M 0.6586628 0.0000000 0.8778808 1.1927859
## Sinkhole.F 1.9027247 0.8778808 0.0000000 2.0927712
## Sinkhole.M 0.7518767 1.1927859 2.0927712 0.0000000
##
## P-values
##
## Marsh.F Marsh.M Sinkhole.F Sinkhole.M
## Marsh.F 1.00000000 0.5101123 0.05707647 0.45212519
## Marsh.M 0.51011230 1.0000000 0.38000838 0.23295324
## Sinkhole.F 0.05707647 0.3800084 1.00000000 0.03636958
## Sinkhole.M 0.45212519 0.2329532 0.03636958 1.00000000
Tests of modularity can also be performed in
geomorph
:
MT <- modularity.test(pupfish$coords,land.gp,CI=FALSE,print.progress = FALSE)
summary(MT)
##
## Call:
## modularity.test(A = pupfish$coords, partition.gp = land.gp, CI = FALSE,
## print.progress = FALSE)
##
##
##
## CR: 0.8771
##
## P-value: 0.003
##
## Effect Size: -4.5427
##
## Based on 1000 random permutations
plot(MT)
One may also evaluate whether the strength of modularity similar in two or more groups. As above this is accomplished using effect sizes based on \(CR\).
coords.gp <- coords.subset(pupfish$coords, group)
modul.tests <- Map(function(x) modularity.test(x, land.gp,print.progress = FALSE), coords.gp)
compare.CR(modul.tests, CR.null = FALSE)
##
## NOTE: more negative effects represent stronger modular signal!
##
##
## Effect sizes
##
## Marsh.F Marsh.M Sinkhole.F Sinkhole.M
## -2.816188 -5.778451 -3.363544 -4.534488
##
## Effect sizes for pairwise differences in CR effect size
##
## Marsh.F Marsh.M Sinkhole.F Sinkhole.M
## Marsh.F 0.000000 2.5466647 1.8303913 1.017783
## Marsh.M 2.546665 0.0000000 0.1496676 1.693290
## Sinkhole.F 1.830391 0.1496676 0.0000000 1.251317
## Sinkhole.M 1.017783 1.6932905 1.2513171 0.000000
##
## P-values
##
## Marsh.F Marsh.M Sinkhole.F Sinkhole.M
## Marsh.F 1.00000000 0.01087579 0.06719144 0.30878104
## Marsh.M 0.01087579 1.00000000 0.88102690 0.09040019
## Sinkhole.F 0.06719144 0.88102690 1.00000000 0.21081881
## Sinkhole.M 0.30878104 0.09040019 0.21081881 1.00000000
One may also wish to determine whether there is support in a given dataset for a 2-module, 3-module, 4-module, or other hypothesis. This is also accomplished via effect sizes of \(CR\).
land.gps3 <- rep('a',56); land.gps3[39:48]<-'b'; land.gps3[c(6:9,28:38)] <- 'c'
#3 module hypothesis (tail now a module)
land.gps4 <- rep('a',56); land.gps4[39:48]<-'b'; land.gps4[c(6:9,28:38)] <- 'c';
land.gps4[c(10,49:56)] <- 'd' #4 module hypothesis (eye now a module)
m3.test <- modularity.test(coords.gp$Marsh.F,land.gps3, iter = 499, print.progress = FALSE)
m4.test <- modularity.test(coords.gp$Marsh.F,land.gps4, iter = 499, print.progress = FALSE)
model.Z <- compare.CR(m3.test,m4.test, CR.null = TRUE)
model.Z
##
## NOTE: more negative effects represent stronger modular signal!
##
##
## Effect sizes
##
## No_Modules m3.test m4.test
## 0.000000 -3.055068 -3.873177
##
## Effect sizes for pairwise differences in CR effect size
##
## No_Modules m3.test m4.test
## No_Modules 0.000000 3.055068 3.873177
## m3.test 3.055068 0.000000 1.030587
## m4.test 3.873177 1.030587 0.000000
##
## P-values
##
## No_Modules m3.test m4.test
## No_Modules 1.0000000000 0.002250096 0.0001074256
## m3.test 0.0022500958 1.000000000 0.3027346572
## m4.test 0.0001074256 0.302734657 1.0000000000