block argument in lm.rrpp / procD.lm / lm.rrpp.ws

[Lucía Alarcón-Ríos]

Dear geomorph community,

I am currently analyzing a morphological dataset (geometric morphometric coordinates) and I have a conceptual question regarding the block argument in lm.rrpp (and procD.lm). I would greatly appreciate your insights.

My sampling design has a hierarchical structure that I need to account for in the analysis:

• My observations (individuals) are naturally clustered into groups (e.g., sampling localities). I have multiple individuals per locality.
• I have two fixed factors: genetic_cluster (2 levels) and colour_morph (3 levels).
• Localities are nested within genetic clusters (i.e., each locality belongs to one of the two clusters).
• All colour morphs occur in both genetic clusters (i.e., the factors are crossed, not nested).

My primary research question is whether colour morphs differ in morphology, while accounting for the underlying genetic structure (cluster) and controlling for allometric effects of size.

If I interpreted it correctly, I need to control for the non-independence of individuals from the same locality (pseudoreplication), as I am not interested in locality as a factor per se. To do this, I have been using the block argument in procD.lm (or lm.rrpp), with the understanding that restricting permutations within localities should preserve the hierarchical structure of the data.

Specifically, my model is:

r

fit06.block <- procD.lm(shape ~ cs * cluster * colour,

                        iter = 9999,

                        RRPP = TRUE,

                        data = mydata,

                        block = locality,

                        SS.type = "II",

                        print.progress = TRUE)

My broadest question is: is this model doing what I think it is doing? That is, is it appropriately accounting for pseudoreplication by restricting residual permutations to occur only within localities?

My doubts arose when I consulted the literature and the package documentation on lm.rrpp.ws. I am aware that this function was designed for other kind of data and questions, but I assumed that lm.rrpp.ws with 'subjects = "locality"' would be conceptually similar to using 'block = locality' in 'lm.rrpp' (well, and SS Type II). However, when I compare the anova tables, I find that the results from lm.rrpp.ws are identical (in terms of Z and significance) to those from lm.rrpp without the block argument. This suggests that 'lm.rrpp.ws' is not applying the same permutation restriction that I expected, which makes me doubt whether I am interpreting the ‘block’ argument correctly.

I have read the relevant literature (Collyer et al., 2015; Collyer & Adams, 2018; Adams & Collyer, 2024) and the package documentation, but it is still not clear for me, probably because my understanding of the underlying statistics is quite limited. I would be very grateful for any clarification you can provide on the conceptual and practical implications of using this model design and structure.

Thank you very much for your time!

Lucía

--

Lucía Alarcón-Ríos, PhD

E-mail: alarconr...@gmail.com
Web: ResearchGate

[Mike Collyer]

Dear Lucía,

Thanks for inquiring about this.  There is a glitch and I will have to resolve it.  I modified the fishy data example for lm.rrpp.ws to mimic your intended analyses and obtained the same kind of baffling results.  My script looked like this:

data(fishy)

suppressWarnings(fit <- lm.rrpp.ws(coords ~  groups * reps,

                                   subjects = "subj", 

                                   data = fishy))

anova(fit)

And sure enough, when I used lm.rrpp with and without blocking, the results were consistent WITHOUT blocking and different WITH blocking.  Not only were they different in terms of Z and P-values, they were different in terms of SS, MS, and Rsq, which should not happen.

However, when I changed the formula above to be this:

suppressWarnings(fit <- lm.rrpp.ws(coords ~  subj + groups * reps,

                                   subjects = "subj", 

                                   data = fishy))

anova(fit)

All ANOVA stats for the non-subject terms were consistent between lm.rrpp.ws and lm.rrpp WITH blocking and SS.type = “II”, which is what one would expect.  Also, the SS, MS, and Rsq were consistent, whether blocking was used, which is reassuring.

So the problem appears to be when a factor is used to establish blocks for permutations but is not part of the formula.  This seems to be one of those cases where an error would be preferred rather than an erroneous result, since the erroneous result slips past our radar.  Because the SS, MS, and Rsq values are not correct when using the block argument and the factor is not in the formula, it suggests a simple problem that cascades into horrible results.  I will investigate it.

In the interim, if you add “locality +” at the beginning of your formula, the ANOVA results for non-locality terms should be trustworthy.  (I confirmed that SS, MS, and Rsq values were consistent for all three models, plus consistent with models not having the subject as a factor, except for the one case where blocking without the term in the formula changed them.). This is inconvenient, as it asks you to ignore the locality term in the ANOVA table, but the non-locality terms should have trustworthy results.

In conclusion, your logic was indeed correct.  The only difference between lm.rrpp.ws and lm.rrpp with a block factor for permutations and SS.type = “II” is that the factor used for subject will have type III SS in the ANOVA table in lm.rrpp.ws; lm.rrpp applies the same SS type to all terms.  Furthermore, lm.rrpp.ws does not allow one to use type I SS.  So you were correct that results should be the consistent with your set-up but you were too kind maybe to think it was your understanding rather than our programming error.

Sorry for this hassle!  I will alert everyone once I resolve this, but it might take a few days.

Best,

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, visit https://groups.google.com/d/msgid/geomorph-r-package/CAHMMSUjhj2%3Ditdh1C6J%3DjpYY2YD-jDgxxWChzOvF3Q5hdA%2BZuw%40mail.gmail.com.