Question on point estimate for CR coefficient in modularity.test

[kieranc...@gmail.com]

While checking the source code for modularity.test, I noted the line CR.obs <- CR(x, gps = gps.obs). It seems that CR.obs represents the CR estimate from the data. If my understanding is correct, the function output should assign CR.obs as the fiinal point estimate for the coefficient. However, in reality, the function used CR.rand[1] as the output point estimate for CR coefficient.

I would like to ask why CR.obs is not directly used as the CR coefficient esitmate? Why assigning the first element of CR.rand, which is random after all. I do not seem to note any line of code that suggests the first element of CR.rand represents the esitmate calculated from real data rather than permuted data.

Thanks.

Kieran 

[Mike Collyer]

Kieran,

The first iteration of any geomorph permutation function is the observed case, which is one random outcome.  (Additionally, one might choose iter = 999 to have 1,000 random outcomes, because the observed case has to count as one.)

You can find this assignment deep in the support code, with the perm.CR.index function (used within the apply.CR function), which creates the iter + 1 sampling frames.  Here is the function code:

function(g, k, iter, seed=NULL){ # g is numeric partition.gp

  if(is.null(seed)) seed = iter else

    if(seed == "random") seed = sample(1:iter,1) else

      if(!is.numeric(seed)) seed = iter

      set.seed(seed)

      p <- length(g)

      ind <- c(list(1:p),(Map(function(x) sample.int(p,p), 1:iter)))

      ind <- Map(function(x) g[x], ind)

      ind <- Map(function(x) as.factor(rep(x,k,each = k, length=p*k)), ind)

      rm(.Random.seed, envir=globalenv())

      ind

}

Notice the bolded part, it is combining a vector of 1:p, for p points with a list of 1:p randomly sampled iter times.  Thus, the first random case is the observed case; i.e., CR.rand[1] is the observed CR.

Hope that helps!

Mike

--
You received this message because you are subscribed to the Google Groups "geomorph R package" group.
To unsubscribe from this group and stop receiving emails from it, send an email to geomorph-r-pack...@googlegroups.com.
To view this discussion on the web, visit https://groups.google.com/d/msgid/geomorph-r-package/60cb624c-6a6c-41f4-b8f8-257e8975caf8n%40googlegroups.com.

[kieranc...@gmail.com]

Thank you for the explanation. Since I am trying to make sure that I really understand CR coefficient, I checked the source code and tried to manually calculate CR coefficient. Hoever, the result I obtained is different from the CR coefficient from modularity.test. I attach my code below and I am puzzled where I was wrong:

library(geomorph)

data(pupfish)

Y.gpa<-gpagen(pupfish$coords, print.progress = FALSE) #GPA-alignment

#landmarks on the body and operculum

land.gps<-rep('a',56); land.gps[39:48]<-'b'

# Method 1: Manual calculation

A <- pupfish$coords

x <- two.d.array(A)

p <- dim(A)[1]

k <- dim(A)[2]

n <- dim(A)[3]

partition.gp <- land.gps

gps <- as.factor(partition.gp)

gps.obs <- as.factor(rep(gps, k, each = k, length = p * k))

g <- gps.obs

gl <- unique(g)

ngps <- length(gl)

gps.combo <- combn(ngps, 2)

Xs <- lapply(1:ngps, function(j) {

            x[, g %in% gl[j]]

    })

ind <- gps.combo[, 1]

a <- ind[1]

b <- ind[2]

S11 <- crossprod(Xs[[a]])

diag(S11) <- 0

S22 <- crossprod(Xs[[b]])

diag(S22) <- 0

S12 <- crossprod(Xs[[a]], Xs[[b]])

sqrt(sum(S12^2)/sqrt(sum(S11^2) * sum(S22^2)))

# Method 2: Use function

MT <- modularity.test(Y.gpa$coords,land.gps,CI=FALSE,iter=99)

My second question is that if I assign each landmark to a separate group, I obtained a CR matrix. I would like to ask if it is acceptable to calculate such point-wise CR coefficient? If there is any cautionary point for such calculation. Thank you.

[Mike Collyer]

It appears that unlike the modularity.test function, you are not obtaining the optimal rotation for the CR statistic.  Because CR is not rotationally invariant and because the orientation of landmark data is arbitrary (for example, whether in GPA data are aligned to the first PC of the landmark configurations), much of the modularity.test function applies to finding the optimal rotation and rotating the data before performing the permutation test.

For the second question, are you asking if it is okay to have, e.g., 10 modules for 10 landmarks?  If that is the question, then the answer is a hard “No!"  Unfortunately, there is no mathematical failure for attempting this, but without within-module covariances (among landmarks), one does not have modules and the CR statistic is meaningless.  If you look at the random values for the statistic in this case, they should all be the same.  We should probably introduce a trap to make sure all modules have at least 2 landmarks.  There might be covariances between x and y coordinates, so unfortunately, the analysis does not blow up.  Regardless, one landmark does not a module make.

Cheers!

Mike

To view this discussion on the web, visit https://groups.google.com/d/msgid/geomorph-r-package/f7afca75-fbed-424a-90ca-d70e94b60120n%40googlegroups.com.

[kieranc...@gmail.com]

Thanks for clarifying a module cannot have only 1 landmark. The 2016 paper that introduced the CR coefficient mentioned that "the CR describes covariation between modules relative to covariation within modules, it is appropriate for evaluating tests of modular structure, rather than simply describing covariation between modules". For my research purpose, I am really interested in pair-wise covariation across the landmarks. I am not particularly interested in between- relative to within-module covariation, as the CR coefficient estimates. Simply describing covariation between modules (in my case: between each pair of landmarks) would satisfy my purpose. I would like to ask how one could "simply describing covariation between" each pair of landmarks?

[Mike Collyer]

Someone else might wish to chime in with a different opinion, but I’ll offer mine since we have been discussing this in this thread.

First, you can obtain a covariance matrix from a gpagen object: $points.VCV is this matrix.

That is a practical answer to your question.  Now for some opinion.  Be very careful with ascribing any biological value to such covariances, as they are not independent of the generalized Procrustes analysis (GPA), which creates covariances based on the spatial distribution and density of landmarks used via the generalized least squares superimposition.  In other words, those covariances have both natural and mechanical sources, with no way of decomposing the two.  One cannot say that “these two landmarks covary more than these other two” and expect that there is a biological implication because of it.  One would have to consider how GPA might have influenced the covariances.

In fact, it is probably a little risky to do the same with modules.  If modules have many points, we might to start to feel more comfortable that stronger covariances within modules rather than between them is evidence of modularity — in spite of GPA — especially if landmarks in different modules are also close together in the coordinate space.  But there is always a risk that meaningful modularity owes some of its effect to GPA and the selection of landmarks for modules.

There are probably better methods to be developed but right now this is the best we have.  You mention “covariation between modules”.  That is integration, a whole different thing, even if the same modules are used.

Cheers!

Mike

To view this discussion on the web, visit https://groups.google.com/d/msgid/geomorph-r-package/1df1f7bf-20da-479c-8766-0a912fe48fd3n%40googlegroups.com.

[kieranc...@gmail.com]

Thank you for the really insightful comments. They are very helpful to me. My research question is on characterizing spatial patterns of covariance of my specimens. I use thousands of spatially dense landmarks to characterize the shape of my specimen. My sample size reaches a thousand. Based on this information, if it is possible to determine how much I could trust findings based on $points.VCV? Or what factors other than sample size and the number of landmarks could influence the validity of the findings?

[lv xiao]

Following Mike's comment that it is inappropriate to assign only 1 landmark per module, I come up with an interesting question.

If the specimen demonstrates object symmetry, is it possible to assign each pair of bilaterally homologous landmarks to the same module (such that the number of modules equals the number of landmarks divided by 2) and then use modularity?

If this can be a solution, which set of procrustes coordinates should be used, the set obtained through gpagen or bilat.symmetry, or both are fine?

To view this discussion on the web, visit https://groups.google.com/d/msgid/geomorph-r-package/070fea9c-a3d3-4fea-a70c-e4d6354eee01n%40googlegroups.com.

[Mike Collyer]

I’m not sure what the goal is, but making sense of a several thousand x several thousand covariance matrix sounds like an interesting challenge.  I recommend using the image function of the Matrix package, e.g.,

library(Matrix)

image(Matrix(mygpa$points.VCV), ...)

This will produce something like a heat map (and can be adjusted by arguments to be specifically a heat map) to show were covariances are strongest or weakest, or most positive or most negative.

Or what factors other than sample size and the number of landmarks could influence the validity of the findings?

I’m not sure what is meant by “findings”.  A covariance matrix is a covariance matrix.  If you mean “validity” insomuch as covariances are biological rather than mechanical, I just don’t think there is anything one can do to make them more biologically valid, other than invent a new superimposition method.

Mike

To view this discussion on the web, visit https://groups.google.com/d/msgid/geomorph-r-package/070fea9c-a3d3-4fea-a70c-e4d6354eee01n%40googlegroups.com.

[Adams, Dean [EEOB]]

As per usual, Mike’s response is spot-on here. 

 

For the current thread, there is a conceptual disconnect between the intention of the CR coefficient and the research question trying to be asked. In short, the two are incommensurate, due in large part to the question being asked, not the analytical tool.

 

Since the 1980s, it has been of interest to evaluate the relative covariation between landmarks. Indeed, this goal seems simple: obtain a valid landmark by coordinate covariance matrix and evaluate the covariance components of interest. However, since the landmark paper (pun intended!) of Rohlf and Slice 1990, it has been known that the superimposition of specimens to standardize position, orientation, and scale, influences the resulting covariance matrix obtained from the aligned coordinates. What this means is that the estimated variance and covariance among landmarks is not accurately represented with this matrix, and so evaluations of it should be treated with extreme caution.

 

The CR coefficient is not a panacea for this problem (see the Discussion of Adams 2016, where there is comment on why this is the case).  

 

The issue is simple. All current analytical methods quantifying modularity parse the signal of some input ‘total’ covariance matrix into that component within modules and that component between modules. To this end the CR ratio does this effectively, and without influence of the number of specimens or variables (Fig 1, Adams 2016). However, the USER must still consider what they are providing the method, and whether it makes sense to do so.

 

For landmark data, we know that many researchers continue to give a GPA-aligned data covariance matrix to some method (CR, RV, etc.) to parse it into modular structure. That may be fine, but one must remember (ala Rohlf and Slice 1990) that some component of the signal is superimposition-induced.  The question then is whether this matters? I simply cannot say: it is data-dependent for now (until we have an alternative superimposition procedure that reduces this effect). But what is certain is that increasing modular structure to the point of having each landmark as its own module will not help matters, and will result in spurious results.

 

In the end one must be a bit more careful in how to define the biological hypotheses being tested, and do so relevant to the methods one wishes to employ. Doing so avoids GIGO.

 

Dean

 

 

Dr. Dean C. Adams

Director of Graduate Education, EEB Program

Professor

Department of Ecology, Evolution, and Organismal Biology

Iowa State University

https://faculty.sites.iastate.edu/dcadams/

phone: 515-294-3834

To view this discussion on the web, visit https://groups.google.com/d/msgid/geomorph-r-package/070fea9c-a3d3-4fea-a70c-e4d6354eee01n%40googlegroups.com.

[Kieran Stringer]

Thank you so much Mike and Dean for help. Based on your replies, I would like to ask is there a way to do some sort of sensitivity analysis to check if research findings and conclusions drawn based on mygpa$points.VCV are robust?

If I have two datasets from two populations with both having a sample size of around 500, and I am not interested in differences in-between. How about obtaining one covariance matrix from each dataset and performing meta-analysis to average each cell in the two covariance matrices.

It seems important that researchers should have some way to gain an understanding of whether their results are overly biased by covariance introduced during GPA.

To view this discussion on the web, visit https://groups.google.com/d/msgid/geomorph-r-package/CH2PR04MB66176506F2A6A922C731E4CDA2549%40CH2PR04MB6617.namprd04.prod.outlook.com.

[lv xiao]

In addition, the question from Lv sounds interesting to learn about. For specimens with object symmetry, if one provides procrustes-aligned symmetric component of the coordinates as the input shape data for modularity.test and assign each pair of bilaterally matching landmarks as a single module, would the resultant CR coefficient valid? There are indeed two points, not one, per module. But the two points are so special that they are symmetric. I cannot tell if there are meaningful within-module variation for CR coefficient to be computed in this case.

To view this discussion on the web, visit https://groups.google.com/d/msgid/geomorph-r-package/CANhDmw_dCnZ6aEmFrbAyYtEb%3Dwi9HVYH3PrMv7owRQGN3Lwb7g%40mail.gmail.com.

[Adams, Dean [EEOB]]

I think that question is inverted. It is not whether the CR makes any sense but rather whether symmetrized data makes sense to use for that application.

 

Dean

 

Dr. Dean C. Adams

Director of Graduate Education, EEB Program

Professor

Department of Ecology, Evolution, and Organismal Biology

Iowa State University

https://faculty.sites.iastate.edu/dcadams/

phone: 515-294-3834

 

To view this discussion on the web, visit https://groups.google.com/d/msgid/geomorph-r-package/CAHCzN%3DQPQL3g5jgg4%2BbSLfiexjcYM9_E0-4N7wi343bVAO-Ycw%40mail.gmail.com.

[kieranc...@gmail.com]

Here is another technical question.

Below is the formula for CR coefficient:

I noted that the line of R code that corresponds to this formula in geomorph is:

sqrt(sum(S12^2)/sqrt(sum(S11^2) * sum(S22^2)))

Since the formula is dealing with trace, I suspect I should see something like sum(diag(S12 %*% S21)) to calculate trace. However, the real code used sum() to calculate the sum of square value of all numbers in the matrix, which I find confusing to me. Similarly, for S11 and S22, since the diagnoal values have been set to 0, squaring the matrix and calculating trace should gave 0 as well, which will make the formula incalculable. I must have missed something but have been puzzled about where I am wrong. Thank you.

[Mike Collyer]

You could look into the method of common principal component analysis:

Flury, B. N. (1984). Common principal components in k groups. Journal of the American Statistical Association, 79(388), 892-898.

There are several ways to compare covariance matrices with CPCA.  

As more of general comment for this entire thread, covariance matrices for shape data are special beasts, not because the covariances are part biological and part mechanical (result of GPA) but because the variables compose a single trait of interest, shape, and not many traits (like landmark locations).  There is a lot of work out there to compare covariance matrices, CPCA being one example, but there is an implicit assumption in most cases that the variables used to estimate covariance matrices have precise meaning, individually, but may covary.  Life history data are good in this way.  Brood size and egg size are two continuous variables with a ratio scale of measurement (0 means 0 and 2 is twice that of 1) that we could analyze independently and expect to negatively covary in some way.  Comparing covariance matrices for data sampled from two populations makes sense, as understanding the correlations among these rather individually important traits might help us infer something about life history strategies.

Doing something like this with shape data is dangerous at best and asinine at worst.  Take for example, Kieran’s example of comparing two covariance matrices for shape data.  Do we lump coordinate data together and perform GPA before splitting them apart again to estimate two covariance matrices or do we perform GPA separately?  If the latter, what if we align configurations to the first specimen in each case?  Arbitrary orientation of the separate specimens will result in different covariance matrices.  Aligning them instead to PC axes will result in different covariance matrices.  Attempting to reconcile these arbitrary choices and find the “biological meaning” in covariances is just impossible, future theoretical advances, notwithstanding.  

But there is no tension here, really.  Kieran states, "results are overly biased by covariance introduced during GPA.”  I assert there is no such thing.  GPA is only a (negatively) biased method insomuch as reducing variation due to size, position, and orientation of landmark configurations is a considered a byproduct rather than a feature.  Sure, call it a bias but it is an intended bias.  Shape is not only a multivariate trait; it is a multidimensional trait, requiring every shape variable in its analysis.  In other words, we might require hundreds or thousands of variables but we have one trait we ascribe as important: shape.  Comparisons of shapes should be the intention when acquiring such data.  Perhaps comparisons of covariance matrices require a different kind of data?

I realize that analyses of modularity and integration feel like covariance matrix comparisons.  I think of them differently.  First, there has to be a natural positive association between landmark density and modularity, at least when modules tend to include proximal neighbor landmarks and exclude distal ones.  Do we really care if adding landmarks makes configurations mode modular?  To my way of thinking, asking if my data shape are modular is boring, but asking how they are modular might be interesting.  Thus, comparing different hypotheses of modularity might be a useful pursuit and I do not think the worry of strong covariances between neighboring landmarks as a result of GPA is a “bias” for which to be concerned.  Rather, examining how the covariances change spatially is a better inference to make.    Integration is greater than expected covariances between hypothetical modules.  Thus, data can be both modular and integrated.  Finding the landscape that bests interprets the pattern of modularity and integration is a good evolutionary goal (especially rather than simply confirming configurations are modular).

There is ample room for methodological improvement, both in terms of measuring and comparing modularity and integration and performing GPA with sliding landmark algorithms.  But if one is asking how can I focus analyses on specific landmarks and their covariances, one has already lost sight perhaps of what shape data are.  (There is currently a discussion about this on Morphmet, as well.)

Cheers!

Mike

To view this discussion on the web, visit https://groups.google.com/d/msgid/geomorph-r-package/CANhDmw_dCnZ6aEmFrbAyYtEb%3Dwi9HVYH3PrMv7owRQGN3Lwb7g%40mail.gmail.com.

[Mike Collyer]

It helps to be a bit more resourceful and do the math before questioning the code.  Here is an example:

> y <- matrix(rnorm(30), 10, 3)

> v <- var(y)

> sum(diag(crossprod(v)))

[1] 7.867005

> sum(v^2)

[1] 7.867005

> library(microbenchmark)

> microbenchmark(sum(diag(crossprod(v))), sum(v^2))

Unit: nanoseconds

                    expr   min      lq     mean  median      uq   max neval cld

 sum(diag(crossprod(v))) 14414 14704.0 15398.54 14871.5 14972.5 56993   100   b

                sum(v^2)   570   588.5   730.09   649.5   758.5  5405   100  a 

So two methods provide the same answer but one is significantly faster and therefore more efficient in thousands of permutations.  Statistics and statistical programming can be different things.

Mike

On May 11, 2021, at 9:28 AM, kieranc...@gmail.com <kieranc...@gmail.com> wrote:

Here is another technical question.

Below is the formula for CR coefficient:

To view this discussion on the web, visit https://groups.google.com/d/msgid/geomorph-r-package/65c42b14-3ae4-4c01-b92a-8d2b6300dae5n%40googlegroups.com.
<Snipaste_2021-05-11_21-16-00.jpg>

[lv xiao]

I revisited this post 4 years ago and noted Dean's statement "Indeed, this goal seems simple: obtain a valid landmark by coordinate covariance matrix and evaluate the covariance components of interest. However, since the landmark paper (pun intended!) of Rohlf and Slice 1990, it has been known that the superimposition of specimens to standardize position, orientation, and scale, influences the resulting covariance matrix obtained from the aligned coordinates." I tried to identify the "landmark paper" by author and year and found one paper entitled "Extensions of the Procrustes method for the optimal superimposition of landmarks" (Syst Zool 39(1):40-59, 1990). However, a brief read of the paper revealed limited explicit discussion of covariance. May I verify this is the "landmark paper"?

Lv

[Adams, Dean [EEOB]]

Yes it is, though the implications are not presented in the manner in which you stated. In this paper it is not discussed as a difference in the estimated covariance matrix, but rather in the variation across landmarks.

See the section "Comparison of Generalized Orthogonal Least-Squares and Resistant-Fit Results" on page 46.  Here it is shown that the input variation from simulation differs from that visualized in the aligned shapes. By definition that affects the estimated PxK covariance matrix.

The discussion section of the paper also has some comments that allude to this (e.g., where they discuss the fact that differences are distributed more evenly over the entire set of landmarks). 

Another early paper that delves into this, and examines the covariance matrices explicity is Walker 2000, which followed up on the observations of Rohlf and Slice. 

Dean

-- 
Dr. Dean C. Adams
Distinguished Professor

Department of Ecology, Evolution, and Organismal Biology
Iowa State University

www.public.iastate.edu/~dcadams/
phone: 515-294-3834

---

From: geomorph-...@googlegroups.com <geomorph-...@googlegroups.com> on behalf of lv xiao <lxi...@gmail.com>
Sent: Tuesday, June 3, 2025 3:47 AM

To view this discussion, visit https://groups.google.com/d/msgid/geomorph-r-package/265f4d58-a821-43e0-bfa2-1a91905eea89n%40googlegroups.com.