Introduction

When one is performing analyses of modularity on sets of traits in different samples, it is often of interest to determine whether those integration values are significantly different from one another. The function compare.CR allows one to statistically compare effect sizes of two or more PLS analyses; in particular, one might wish to compare modularity between two or more samples, each measuring modularity between separate sets of traits. Alternatively, this method can also be used to compare the degree of modular signal among alternative hypotheses of modularity for the same dataset.

This analysis calculates effect sizes as standard deviates, z, and performs two-sample z-tests, using the pooled standard error from the sampling distributions of the PLS analyses (Adams & Collyer, 2019).

The input for this function must be of class “pls,” that is, an object resulting from the functions modularity.test, or phylo.modularity. Any number of objects can be input.

compare.CR()
  • \(...\): saved analyses of class CR
  • \(CR.null\): A logical (TRUE/FALSE) value to indicate whether a Null CR model (no modularity) should also be included in analysis. When comparing alternative hypotheses of modularity, this should be set to TRUE.
  • \(two.tailed\): A logical value to indicate whether a two-tailed test (typical and default) should be performed.



Comparing Modular Signal Across Datasets

By way of example, we present an analysis of modularity between the body and operculum of pupfish. These data are included with geomorph by default, as a single object. The data have been separated by population (marsh vs. sinkhole) and sex for the purposes of this example. Please see the tutorials on Data Manipulation and R Data Basic for information on how to accomplish this.

Here, the object ‘group,’ refers to the factor used to partition the pupfish data.

levels(group)
## [1] "Marsh.F"    "Marsh.M"    "Sinkhole.F" "Sinkhole.M"



Now, we run our modularity analyses with the help of the map function. This allows use to run the function on each of the partitioned sets of coordinates.

modul.tests <- Map(function(x) modularity.test(x, land.gps,iter=999, print.progress = FALSE), coords.gp)

This returns a list object containing the results of our four integration tests. Remember that objects within a list can be accessed using the $ operator.

Finally, we perform the statistical comparison of our results.

group.Z <- compare.CR(modul.tests, CR.null = FALSE)
summary(group.Z)
## 
##  NOTE: more negative effects represent stronger modular signal! 
## 
## 
## Effect sizes
## 
##    Marsh.F    Marsh.M Sinkhole.F Sinkhole.M 
## -0.7265087 -4.1650804 -2.4916997 -2.1320218 
## 
## Effect sizes for pairwise differences in CR effect size
## 
##             Marsh.F    Marsh.M Sinkhole.F Sinkhole.M
## Marsh.F    0.000000 2.45689940 1.74666288  1.0170828
## Marsh.M    2.456899 0.00000000 0.07932252  1.4191065
## Sinkhole.F 1.746663 0.07932252 0.00000000  0.9811416
## Sinkhole.M 1.017083 1.41910651 0.98114159  0.0000000
## 
## P-values
## 
##               Marsh.F    Marsh.M Sinkhole.F Sinkhole.M
## Marsh.F    1.00000000 0.01401419 0.08069583  0.3091141
## Marsh.M    0.01401419 1.00000000 0.93677609  0.1558680
## Sinkhole.F 0.08069583 0.93677609 1.00000000  0.3265229
## Sinkhole.M 0.30911405 0.15586797 0.32652292  1.0000000

Summarizing these results, returns two tables of pairwise z-tests and p-values, as well as the effect sizes for each of the modularity analyses.

Compare Alternative Modular Hypotheses

Finally, we illustrate here how one might compare alternative hypotheses of modularity. For this example, we have modularity analyses on two separate partitions of the same data; one that hypotheses three modules (land.gps3), and one that hypotheses four modules (land.gps4). The process is similar to above. Note that we are running these analyses on female individuals only.

First, we run out modularity analyses:

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)

Then compare them. Note that, since we are comparing alternate hypotheses, the CR.null argument is set to TRUE.

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

summary(model.Z)
## 
##  NOTE: more negative effects represent stronger modular signal! 
## 
## 
## Effect sizes
## 
##          No_Modules modul.tests$Marsh.F             m3.test             m4.test 
##           0.0000000          -0.7265087          -2.6310761          -3.8848807 
## 
## Effect sizes for pairwise differences in CR effect size
## 
##                     No_Modules modul.tests$Marsh.F   m3.test   m4.test
## No_Modules           0.0000000           0.7265087 2.6310761 3.8848807
## modul.tests$Marsh.F  0.7265087           0.0000000 1.2695040 2.0305221
## m3.test              2.6310761           1.2695040 0.0000000 0.7533947
## m4.test              3.8848807           2.0305221 0.7533947 0.0000000
## 
## P-values
## 
##                       No_Modules modul.tests$Marsh.F     m3.test      m4.test
## No_Modules          1.0000000000           0.4675270 0.008511496 0.0001023802
## modul.tests$Marsh.F 0.4675269581           1.0000000 0.204261351 0.0423034965
## m3.test             0.0085114962           0.2042614 1.000000000 0.4512127816
## m4.test             0.0001023802           0.0423035 0.451212782 1.0000000000

The result of a summary in this case, returns tables of pairwise effect sizes and P-values. This also includes the null hypothesis of no modularity.

This function returns an object of class “compare.CR”, which is a list containing the following:

compare.CR output
  • \(sample.z\): A vector of effect sizes for each sample.
  • \(sample.r.sd\): A vector of standard deviations for each sampling distribution (following Box-Cox transformation).
  • \(pairwise.z\): A matrix of pairwise, two-sample z scores between all pairs of effect sizes.
  • \(pairwise.p\): A matrix of corresponding P-values.