Multiple comparisions between sexes within years ANOVA or t-tests? Tukey posthoc? Bonferroni?

Question

I am conducting a genetic analysis comparing Fst and relatedness coefficient values between male and female animals across years from 2013-2019. Each year is divided into quarters, so I have a total of 40 quarter-year categories. Each year has a varying number of individuals, and the number of males and females are not equal in each year. I want to see if values for the sexes differ within and between each quarter year. I am using R so I did a linear model (Fst ~ Sex*Year) and then did ANOVA on the model with a Tukey posthoc test to see the comparisions between the sexes within each quarter year.

Fst.test = lm(Fst ~ Sex*Year, data = Fst)
Fst.aov = aov(Fst.test)
summary(Fst.aov)
Fst.posthoc = TukeyHSD(Fst.aov)

Many of my posthoc comparisions of the sexes within each quarter-year were not significantly different, even though the values looks quite different. So I chose a few at random and did a standard t-test, and they were signficiant that way.

What stats should I be doing here? Is ANOVA and Tukey posthoc tests the best way or should I be doing multiple t-tests for each quarter with some kind of correction? If so, which correction, I remember doing a Bonferroni years ago...

Answer 1

First, as I understand it, the genetic fixation index $F_{st}$ is limited to values between 0 and 1, so a simple linear regression model might not be applicable. I suspect that similar work in this field might have used other forms of regression that are designed to handle outcomes over this limited range, like beta regression. The Wikipedia page linked above has links to some R implementations of software for $F_{st}$, which might provide subject-matter-specific guidance.

Whichever regression model you end up using, however, you are asking a lot by trying to do 40 multiple comparisons at once. Any correction for the multiple comparisons problem is, appropriately, going to give you higher p-value estimates than you would get from single t-tests. To make this concrete, the (unnecessarily conservative) Bonferroni correction will require $p<0.00125=0.05/40$ for a single comparison to get a family-wise error rate (Type I error) of 0.05 with 40 comparisons. If acceptable in your field, you might correct for the "false discovery rate" instead (e.g., Benjamini-Hochberg; see the R p.adjust() function), accepting that a certain fraction of your "true" differences are actually going to be false.

Do you need all 40 comparisons? You might be better off treating year as a continuous predictor fit flexibly with a spline (perhaps including a seasonal variable), which would require many fewer coefficient estimates and fewer comparisons. The emmeans package provides tools for post-modeling analysis and multiple-comparison correction for a much broader class of models and choices of comparisons than the basic TukeyHSD() function can handle.

There might also be some correlations over time that aren't adequately handled by the independence assumption inherent in your linear model; presumably the $F_{st}$ at one observation time is correlated via inheritance to prior $F_{st}$ values. This is outside my particular expertise; consult with a practitioner experienced in field genetics if you don't get more help here.