Null model in pairwise comparisons

[Lucía Alarcón Ríos]

Dear geomorph community,

I would like to ask a question regarding the choice of the null model in function pairwise in RRPP and its implications in results interpretation.

I have a factorial model: (shape ~ f1*f2) with a significant interaction.

I guess that the accurate null model to explore the interaction term would be (shape ~ f1 + f2). However, I am interested in differences among groups (i.e, groups = f1 x f2), thus, I wonder whether using a null model including only the intercept (i.e., shape ~1) would be appropriate.

Indeed, I am aware that the latter null model implies that pairwise comparisons would include general shape variation due to f1 and f2 and their interaction, but I do not fully understand the implications in results interpretation.

Thank you in advance,

Best

Lucía

--

Lucía Alarcón-Ríos

Margarita Salas Postdoctoral Researcher

1- Universidad de Oviedo

2- Visiting Researcher at Animal Ecology Group (University of Vigo) / Investigadora visitante en Grupo de Ecoloxía Animal (Universidad de Vigo)

E-mail: alarco...@uniovi.es  / lucia.ala...@uvigo.es / alarconr...@gmail.com
Web: ResearchGate

[Ian Dworkin]

Hi Lucía,

A few years ago at the beginning of the pandemic we had a virtual morphometrics meeting, and I did this introduction (sort of a deep dive) into using this function,

Here is a link to the code and notes to that session.

https://github.com/DworkinLab/VirtualMorphoMeetup#using-pairwise-in-geomorph-id

I think we recorded the video of the zoom session and I will see if I can find that link as well.

Cheers

Ian

--
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/CAHMMSUh7m0hdXhnYdLRYYZCRkLAoDd356xwiohwYNWLZpwfedw%40mail.gmail.com.

--

Ian Dworkin
Department of Biology
McMaster University

Office phone 905 525 9140 ext. 21775

Lab phone 905 525 9140 ext. 20076
dwo...@mcmaster.ca

dworkinlab.github.io

[Lucía Alarcón Ríos]

Hi Ian,

Thank you for your response and the link to your tutorial. I found it really helpful.

According to the tutorial, using a null model that only includes the intercept would be useful to " to ask whether the means between pops (f1) are different as well as sex (f2)". But, a few lines below you pointed that " Perhaps the appropriate question given the ANOVA results and the "PC" plot is whether there is evidence of "significant" SShD in either population AFTER accounting for additive effects of Sex and Pop."

Precisely, that was my main concern. Considering that I have a significant interaction, I would like to know whether comparing the overall effects (f1, f2, and f1:f2) by using only a null model that includes the intercept would be theoretically reasonable.

Thanks!

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

[Mike Collyer]

Dear Lucía

Be careful not to conflate null models for random data generation and null models for testing in your logic. You can use a model with only an intercept to generate null distributions but in ANOVA, every term compares two models’ sums of squares and cross-products (SSCP)— a reduced model and a full model — with the difference between them being the inclusion of the term that is tested. One can use either type I, type II, or type III SSCP, but for a two-factor factorial model and the typical way ANOVA is performed, they would all be the same. The reduced model would contains f1 and f2, the full model both effects plus their interaction. Thus, in this kind of ANOVA there is no way to state it differently. The test of the interaction accounts for the main effects (additive if using type I SSCP but it would not matter if it was not additive for this one effect).   But there are multiple ways to perform ANOVA. 

If you want to test all groups, fit two models:

fit.null <- procD.lm(shape ~ 1, …)

fit.full<- procD.lm(shape ~ f1 * f2, …)

And test this with

anova(fit.null, fit.full, …)

groups = interaction(f1, f2)

pairwise(fit.full, fit.null = fit.null, groups = groups)

The anova and pairwise tests here are commensurate and use an intercept model as a null. The anova output provides p-values and z-scores assuming a null model for only an intercept. The important step for ANOVA is to have both models as arguments to override what the null model should be. 

Hope that helps,

Mike

Sent from my iPhone

On Oct 11, 2022, at 2:25 AM, Lucía Alarcón Ríos <alarconr...@gmail.com> wrote:



To view this discussion on the web, visit https://groups.google.com/d/msgid/geomorph-r-package/CAHMMSUittTPF5PuszHt5bCxHQXORB-zezSs_Wk3o5oWMYqdC%2BA%40mail.gmail.com.

[Ian Dworkin]

Dear Lucia,

As Mike mentioned, if you are interested in examining differences in all combinations of groups from your model (shape ~ f1 + f2 + f1:f2) then the appropriate null model is (shape  ~ 1). 

In the tutorial we were largely focused on an example where I framed it as us only interested in sexual shape differences. So we wanted a reduced model (which I generally prefer for terminology instead of "null") where we accounted for one factor (in that case population). So we were effectively asking questions like "what are the sex differences after having accounted for the effects of population?"

One thing to be aware of though is that if you do not have 1-2 specific planned comparisons to test your predictions among all the combinations of f1 and f2, and you are really "examining" all combinations of f1 by f2, then you may have a severe multiple testing problem (depending on total number of groups) that you want to account for. pairwise is "agnostic" to such issues and expects you to make the appropriate judgement calls of how to do so.  So something like a Bonferroni correction (if you do not have a very large number of comparisons) among other approaches may be appropriate.

I will say that I try to steer people away from all by all comparisons, and instead to focus on a priori tests of specific predictions, unless there are very good reasons not to.

I apologize that I was unable to find the zoom recording. I think it may have been during my first times on zoom (hard to imagine now) and I did not record that session or something? Not sure, but I am unable to find it. I think that those of us chatting discussed some of these very issues in more detail.

All the best

Ian

To view this discussion on the web, visit https://groups.google.com/d/msgid/geomorph-r-package/B04E38A9-AA3B-4BA9-AE91-B5BC48F6A040%40gmail.com.

[Lucía Alarcón Ríos]

Dear Mike and Ian, thank you very much for your answers.

I would like to clarify a few things from your comments, just to be sure that I have understood the rationale of your explanations. Sorry if there is something inaccurate or very repetitive because of my lack of statistical knowledge.

1) The anova test suggested by Mike:

anova(fit.null, fit.full, …)

I have always thought that this was used to compare models (I'm missing something for sure), how does it fit with the specification of the null model you want to use?. I mean, does it imply that the test of every term in the full model 'fit.full' (i.e., f1, f2 and f1:f2) is performed accounting for a reduced model that only contains the intercept and, after this, both models (intercept only and full) are compared between them?

2) Ian, I am curious precisely about the example you mentioned from the tutorial where you tested differences among sexes after considering population effects. If I'm right, because you have a significant interaction you cannot just 'remove' population variation. And now, I think I'm going to answer myself, that's why you add a note pointing out that "We should really (given our hypothesis) only look at the the M-F differences within each population"?

3) Finally, regarding multiple comparisons. I have 6 groups (15 pairwise comparison), and because of the issues that you mentioned, I decided to apply sequential Bonferroni corrections. However, do you consider that subsetting data for more directed pairwise comparisons is a better approach? Actually, I just tried this approach and I get the same results either with all comparisons+sequential Bonferroni and subsetting,

Again thank you for your fruitful comments and explanations,

Best,

Lucía

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

[Mike Collyer]

Dear Lucía,

1) The anova test suggested by Mike:

anova(fit.null, fit.full, …) 

I have always thought that this was used to compare models (I'm missing something for sure), how does it fit with the specification of the null model you want to use?. I mean, does it imply that the test of every term in the full model 'fit.full' (i.e., f1, f2 and f1:f2) is performed accounting for a reduced model that only contains the intercept and, after this, both models (intercept only and full) are compared between them? 

Yes, it is to compare models, and any ANOVA analysis compares models.  What one might normally do, for example,

anova(fit.full)

is a series of model comparisons, specifically the SS of competing models.  When only one model fit is supplied, the function attempts to find the series of model comparisons to use; if two or model fits are used, the first one in is considered the null/reduced model for the others. Using your example, shape ~ f1 + f2 + f1:f2, you might get ANOVA tables with effects estimated for f1, f2, and f1:f2.  If SS = type I, then here is the series of model comparisons:

effect reduced full

f1 intercept intercept + f1

f2 intercept + f1 intercept + f1 + f2

f1:f2  intercept + f1 + f2 intercept + f1 + f2 + f1:f2

If SS = type 2

effect reduced full

f1 intercept + f2 intercept + f1 + f2

f2 intercept + f1 intercept + f1 + f2

f1:f2  intercept + f1 + f2 intercept + f1 + f2 + f1:f2

This is to say that every line in an ANOVA table is basically an application of anova(reduced model, full model), based on a strategy dictated by the type of SS.  The only exception is that the residual SS in the F-statistic uses the full-model residuals.

By using anova(fit.null, fit.full, …), you are directing the analysis to perform a comparison as such:

effect reduced full

f1 + f2 + f1:f2 intercept intercept + f1 + f2 + f1:f2

Here is something else you could do.  You could make a new variable,

groups = interaction(f1, f2)

and a new model fit

fit.g <- procD.lm(shape ~ groups, …)

and test this with ANOVA:

anova(fit.g)

and you will get the same exact result as anova(fit.null, fit.full, …) or as anova(fit.null, fit.g)  (The effect of group is the same as the combined effects of f1, f2, and f1:f2.)  A factorial model can be made into a single-factor model (but as Ian noted, perhaps with some concerns).  With a single-factor model, you only have one model comparison to make, so strategies converge.

I think this is the same thing you were asking and I hope my response was not overly pedantic, but I hope this might help anyone reading this thread who is trying to wrap their head around these subtle but similar options.

Cheers!

Mike

[Ian Dworkin]

Lucía,

The anova test in RRPP/geomorph always has some form of reduced model(s). You can get a sense (at least with respect to Type I sum of squares) what it defaults to with the reveal.model.design(). I usually suggest specifying the reduced model yourself so you know exactly what you are comparing to. While RRPP implements this using a permutation test, it is essentially doing a classic ANOVA test between the models (see slides 18-26 or so in the R markdown ones you knit it into slides... I go over those points there). If you specify other types of SS, then there will be multiple reduced models (depending on what model is being compared) to assess particular parameters in the model. ... Mike just responded to this, so I will stop this part....

For your second question point "If I'm right, because you have a significant interaction you cannot just 'remove' population variation.". Well not really in a general sense...

Let me give you a completely made up example with the pupfish data set. Let's say that a priori we were aware that shape was a sexually selected trait in this species. Furthermore, one of the populations had a predator species present and experienced high predation rates.  We had some evidence (from other data, independent of what we are currently examining) that variation in body shape influenced performance (swimming speed or something like that) and ultimately the ability to get away from the predators (i.e. shape is a phenotypic target of viability selection in the population with the predators). So we predict that we would observe changes in the amount of sexual shape dimorphism, and that shape as a sexual signal may be reduced (potentially leading to a reduction in sexual shape dimorphism), but only in the one population.

So here we a priori know that there should be sexual shape dimorphism (from previous work) and from the previous study, likely population differences in shape (due to predation). However what we want to know, and have collected this data to test is whether sexual shape dimorphism is the degree to which sexual shape dimorphism differs between these populations. So in this case in the context of the ANOVA what we are interested in is the interaction term. Similarly in the pairwise contrasts we are primarily interested in whether the magnitude of the sexual shape dimorphism has changed (and changed in a manner consistent with our expectation based on predation). In this case I would use a full model of  (shape ~ sex + pop + sex:pop) and a reduced model of (shape ~ sex + pop) so that in my use of pairwise I am specifically addressing the hypothesis of change (after having accounted for the additive effects of both sex and population on shape) in sex. In other words whether the changes in shape in the data (in particular for the population I am interested in) are different (in this case less than) what we would expect based on the overall differences of population and sex. There are actually other ways of addressing this hypothesis a bit more cleanly, but in terms of pairwise, this would be a reasonable approach.  Hopefully this rambling example makes sense.

For your third point, rarely do I recommend doing analyses on subsets. Generally fit the full model on the complete data set and set up your specific contrasts if you are able to. While it is too late (since you have already done the analysis), in a future study with the 15 possible comparisons, if you a priori write down the specific comparisons you are interested in, then you can either (using lm or emmeans) set up the specific contrasts, or with pairwise, only examine the ones you are interested (and still use Sequential bonferroni or other approaches). The reason to use all of the data and the full model (instead of subsets) usually comes down to more appropriate ways of estimating uncertainty in your model estimates. If you have a perfectly balanced design (and all factors) then it won't matter, but for most datasets, imbalanced data generates a form of co-linearity, meaning that the various parameters associated with each factor in your model can not be estimated completely independently of one another.

Hope that all made sense. I am replying quickly between meetings, so apologies for any typos.

Ian

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

[Lucía Alarcón Ríos]

Dear Mike,

Thank you very much for your insightful explanation.

Now I can say that those basic and essential concepts are much more clear for me.

I strongly recommend reading this to everyone.

Thank you!

Lucía

--
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/AF928443-AF63-467C-AB6E-D530F52A123D%40gmail.com.

[Lucía Alarcón Ríos]

Dear Ian,

Your example is now much  more clear. I think we were talking about different examples from the tutorial (I was referring to the one with a reduced model 'shape~1+pop' - in Line 391). But your detailed explanation was really helpful to confirm my thoughts and that I had understood it properly.

Thank you also for your suggestions about multiple pairwise comparisons. I will try to implement it into my analysis, as I found it really appropriate and it is something that I had thought about a lot, because it is a problem that always arises when working with multiple comparisons.

Thank you for your help,

Lucía

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

[Hardeep Kaur]

Dear Ian, Mike, and Lucia,

Thank you so much for this fruitful thread.

Kind regards,

Hardeep Kaur

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