I have an unbalanced design with nested factors. In particular, the nested factor is obesity, with levels: of obese and nonobese. The other factor is group: women with syndrome, control women, and control men. The response let's say is a metabolite.
What I have done, is to set a linear model, with the following contrasts, due to the nature of the problem: the coefficients of the variates sum up to zero.
Then, in R one can achieve this with:
mdl <- lm(Y ~ obesity*group,
contrasts=list(obesity=contr. sum,
group=contr.sum)
car::Anova(mdl, type=3)
Then the output is
Anova Table (Type III tests)
Response: Y
Sum Sq Df F value Pr(>F)
(Intercept) 549.02 1 1117.2302 < 2e-16 ***
obesity 0.19 1 0.3927 0.53446
group 2.09 2 2.1311 0.13199
obesity:group 3.21 2 3.2632 0.04866 *
Residuals 19.66 40
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
I tested for the assumptions of normality of the residuals and homogeneity, the model accomplishes this.
Now, as you can see, the interaction is significant, however, running TukeyHSD test, is not well suited because, from the help in R, it says:
" This is a generic function: the description here applies to the method for fits of class "aov" "
Thus, aov functions are the wrapper of anova(lm), but it does not take into consideration the contrasts, and it only does the sequential sum of squares
So I think, that the problem could be resolved with "marginal means"
The problem: I have no idea how to insert the formula with the error term.
I tested:
mdl2 <- aov_car(Y ~ obesity*group +
Error(patients)
Where patients are the subjects. I know that the error term, is the random mixed effect, (I think..)
Then I can run
emmeans(mdl2)
with the same output,
But I cannot figure out which is the significant group of interaction that the above results say.
So the question is why is significant the ANOVA but not the interaction ? Should I run pairwise t test with interaction ? I mean, interaction(obesity,group) and then run pairwise t.test ?
I have seen several posts about this issue, in particular, I think that the question is already answered on this package: https://cran.r-project.org/web/packages/emmeans/index.html
More than the data itself, it's a question of statistics