I can't work out if I should be using 2 or 3 factors here.

  • I have two groups: CON and TEST
  • I gave Drug A or Drug B (there is no "None" group...)
  • I take samples before treatment and after treatment (before and after) (so technically the "before" = "None")


  • If I do 2 x 2 (Group x Drug) analysis, then I can investigate the effects of group, but I won't know which drug is having an effect if there is an effect of drug.

  • If I do a 2 x 2 x 2 (Group x Drug x Time) analysis, I worry that drug and time aren't independent.

I want to know if there is an effect of group, an effect of either drug, and specifically an effect of Drug B and if there are any interactions. Would it be more sensible to do a 2 x 2 and then a post-hoc test in this case?

In practice I am using a mixed effects model with "Sample ID" as a random effect to account for the repeated sampling (before and after).

 fit <- nlme::lme(Measurement ~ Group*Drug, random = ~1|SampleID, method="REML", data = data.df)

Edit: I've created an example df to see the groups

   Drug  Group Time   SampleID Measurement
   <fct> <fct> <fct>     <int>       <int>
 1 A     Con   Before        1           3
 2 A     Con   After         1           5
 3 A     Con   Before        2           5
 4 A     Con   After         2           1
 5 A     Test  Before        3           5
 6 A     Test  After         3           4
 7 A     Test  Before        4           1
 8 A     Test  After         4           4
 9 B     Con   Before        5           2
10 B     Con   After         5           2
11 B     Con   Before        6           5
12 B     Con   After         6           3
13 B     Test  Before        7           3
14 B     Test  After         7           3
15 B     Test  Before        8           4
16 B     Test  After         8           1

2 Answers 2


You should fit the model consistent with your design. In this is case, you have 3 factors, so you need a 3-way ANOVA. It doesn't matter if you don't have the power to detect the 3-way interaction; the correct model should include all interactions that might be present. You can use planned contrasts to ask the specific questions you want to answer, especially if your goal is to compare specific cells or combinations of cells.

You don't need a multilevel model for this; you can use a regular linear model and use cluster-robust standard errors or cluster-bootstrap standard errors. This design is also specifically suited for an analysis called Mixed ANOVA. Assuming you go with a random effects model, it would look like the following:

fit <- lmerTest::lmer(Measurement ~ Group * Drug * Time + (1|SampleID), data = data)

You can then follow it up by using the tools in the marginaleffects package to probe specific differences. For example, to examine the effect of group (marginalizing over other variables), you can run

avg_comparisons(fit, variables = "Group")

To examine the effect of time within each drug (marginalizing over group), you can run

avg_comparisons(fit, variables = "Time", by = "Drug")

To examine the effect of Time within each level of group and drug, you can run

avg_comparisons(fit, variables = "Time", by = c("Group", "Drug"))

To see whether the effect of Time differs between Drug A and Drug B within the TEST group, you can run

avg_comparisons(fit, variables = "Time", by = c("Group", "Drug"),
                hypothesis = "b4 - b3 = 0")

Essentially, any quantity you want to compute that is a function of the group means can be computed using this framework. Note that the questions you want to answer are separate from the model you use. You can run all the same functions regardless of how you parameterize the model. You should not parameterize your model to have interpretable coefficients; rather, you should fit the model that best characterizes your data and use planned contrasts like these above to answer the questions you have about your quantities of interest.

  • $\begingroup$ Hmm perhaps the wording in my answer was not great. I agree that you should fit the model to what your specified theory should be, I was more warning that interactions can be data hungry. Your advice here is a more complete answer in any case (+1) $\endgroup$ Feb 24 at 0:59

It looks to me that you have a straightforward two-way factorial analysis that can be fit to a mixed model quite easily (using lme4 syntax below):

fit <- lmer(y ~ Factor_A * Factor_B + (1|ID), data = data)

This syntax estimates the differences in mean values of $y$ for the independent effects of both factors, their interaction, and the random intercept of the repeated measures (IDs of the people involved). Note that this just estimates two clusters for the random effect, so it is entirely plausible that there may not be enough variation between time points (which may lead to a failure to converge). If this happens, one alternative is to instead model it as a fixed effect to account for the autocorrelation in some way (though doing so as an interaction will result in a three-way interaction that may not be powerful enough to detect effects if your sample is small).

  • $\begingroup$ Thank you for your quick answer! I've added extra details in the question to ask this: in a 2x2 design, wouldn't that pick up any "effects" of Drug in the Before groups? I.e. the "before" groups would contribute the any effects of Drug, even though they haven't yet received either drug. Hope that makes sense? Also, if I pass this to emmeans it wouldn't give me useful comparisons in regard to the before/ after groups $\endgroup$
    – DS14
    Feb 23 at 13:55
  • $\begingroup$ If you include the "before" groups, then yes, you would include the baseline effects before treatment. That is not a bad thing. You ideally want something to compare against, to see if the effect is indeed null or substantial. But this again depends on if you model it as fixed or random. This part of the model would essentially be missing from pairwise tests if it is not modeled as fixed. Whereas adding it in as fixed will directly determine the conditional mean differences at both time points. $\endgroup$ Feb 23 at 14:12
  • $\begingroup$ Okay. So in that case I would do: fit <- nlme::lme(Measurement ~ Group*Drug + Treatment, random = ~1|SampleID, method="REML", data = data.df) or should it be an interaction term (Group*Drug*Treatment)? $\endgroup$
    – DS14
    Feb 23 at 14:17
  • $\begingroup$ You would essentially use the first model if you are not interested in time effects from a deterministic perspective and the second model if you are determined to see if drug means by group differ with respect to the time of treatment. If one only wants to include the temporal aspect independent of the other effects, then a main effect of treatment without an interaction can also be specified. $\endgroup$ Feb 23 at 14:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.