# Compute ANOVA F if you only know post-hoc results

It's not uncommon to read a scientific article in which the authors perform an ANOVA followed by post-hoc t-tests, but don't report the results of the ANOVA at all, and instead skip straight into the results of the post-hoc t-tests. Assuming all post-hocs are reported, is there any way that a reader could back-calculate the results of the ANOVA? For simplicity, let's say it's a 2-way ANOVA.

Edit here is an example:

Load data

library(palmerpenguins)
library(emmeans)

data(package = 'palmerpenguins')
dat = penguins


Categorical x continuous interaction

model1 = lm(body_mass_g  ~ species*flipper_length_mm ,data=dat)
# anova
anova(model1)
Analysis of Variance Table

Response: body_mass_g
Df    Sum Sq  Mean Sq F value    Pr(>F)
species                     2 146864214 73432107 534.660 < 2.2e-16 ***
flipper_length_mm           1  24776495 24776495 180.398 < 2.2e-16 ***
species:flipper_length_mm   2   1519564   759782   5.532  0.004327 **
Residuals                 336  46147424   137344
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# post-hoc tests
res1 = emtrends(model1, pairwise ~ species, var = "flipper_length_mm")
res1$contrasts contrast estimate SE df t.ratio p.value Adelie - Chinstrap -1.74 7.86 336 -0.222 0.9733 Adelie - Gentoo -21.79 6.94 336 -3.139 0.0052 Chinstrap - Gentoo -20.05 8.19 336 -2.448 0.0394 P value adjustment: tukey method for comparing a family of 3 estimates  Categorical x categorical interaction model2 = lm(body_mass_g ~ species*sex,data=dat) # anova anova(model2) Analysis of Variance Table Response: body_mass_g Df Sum Sq Mean Sq F value Pr(>F) species 2 145190219 72595110 758.358 < 2.2e-16 *** sex 1 37090262 37090262 387.460 < 2.2e-16 *** species:sex 2 1676557 838278 8.757 0.0001973 *** Residuals 327 31302628 95727 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 # post-hoc tests res2 = emmeans(object = model2, pairwise ~ species|sex) res2$contrasts
sex = female:
contrast           estimate   SE  df t.ratio p.value
Adelie - Chinstrap     -158 64.2 327  -2.465  0.0377
Adelie - Gentoo       -1311 54.4 327 -24.088  <.0001
Chinstrap - Gentoo    -1153 66.8 327 -17.246  <.0001

sex = male:
contrast           estimate   SE  df t.ratio p.value
Adelie - Chinstrap      105 64.2 327   1.627  0.2357
Adelie - Gentoo       -1441 53.7 327 -26.855  <.0001
Chinstrap - Gentoo    -1546 66.2 327 -23.345  <.0001

P value adjustment: tukey method for comparing a family of 3 estimates


In either case, is it possible to compute any of the ANOVA results, using only the post-hoc test results?

• Please post an example dataset.csv, some (simple!) ANOVA results, and t-test results. If you also include the code you used, that would be helpful to folks responding.
– J_H
Jan 8 at 2:46
• I've added some examples. Jan 8 at 16:28

## 1 Answer

In principle, this may be possible if you have a completely randomized design. When there are blocks, pairing, or mixed effects, it gets harder or impossible.

Simplest case is with a table of cell means and their standard errors:

> emmeans(model2, ~sex*species)
sex    species   emmean   SE  df lower.CL upper.CL
female Adelie      3369 36.2 327     3298     3440
male   Adelie      4043 36.2 327     3972     4115
female Chinstrap   3527 53.1 327     3423     3632
male   Chinstrap   3939 53.1 327     3835     4043
female Gentoo      4680 40.6 327     4600     4760
male   Gentoo      5485 39.6 327     5407     5563

Confidence level used: 0.95


Then you can use the emmobj function to construct an emmGrid object:

> eobj = emmobj(c(3369,4043,3527,3939,4680,5485), V=diag(c(36.2,36.2,53.1,53.1,40.6,39.6)^2),
+     levels = list(sex=c("f","m"), species=c("a","c","g")), df = 327)

> eobj # shows same table as before
sex species estimate   SE  df lower.CL upper.CL
f   a           3369 36.2 327     3298     3440
m   a           4043 36.2 327     3972     4114
f   c           3527 53.1 327     3423     3631
m   c           3939 53.1 327     3835     4043
f   g           4680 40.6 327     4600     4760
m   g           5485 39.6 327     5407     5563

Confidence level used: 0.95

> joint_tests(eobj)
model term  df1 df2 F.ratio p.value
sex           1 327 311.580  <.0001
species       2 327 747.804  <.0001
sex:species   2 327   8.735  0.0002


This is comparable to the Type III anova for the model:

> joint_tests(model2)
model term  df1 df2 F.ratio p.value
species       2 327 746.924  <.0001
sex           1 327 311.838  <.0001
species:sex   2 327   8.757  0.0002


But the anova you show is the Type I (sequential) anova, and you can't get that from eobj.

Now, with only the pairwise comparisons shown in the illustration, you do't have any information about the sex effects or sex:species effects. If you had also the results for pairwise comparisons of sex for each species, I think it may be possible to cobble a table of means and SEs, though with a different grand mean, by arbitrarily setting one of the cell means to zero and working from there. But it gets messy because each SE you see is $$\sqrt{1/2}$$ times the square root of the sum of two cell variances. In principle, it's possible to work through this but it won't be fun.

It is not possible to do the anova for model1 where we have comparisons of slopes. Again it may be possible if we know the slope estimates themselves, and some kind of information about the interaction. But it would be very ugly.

Again, I emphasize these ideas won't work for mixed models or others where cell means are in any way correlated.