Correct approach to comparing sub-group means across multiple treatment conditions

Question

Let's say I have an experiment where I randomly pair people up with either a male or a female player, and they play a game together. I then measure some dependent variable -- call it outcome . From the perspective of participants, there are four possible subgroups -- Male-Male (MM), Male-Female (MF), Female-Male (FM), and Female-Female (FF). I have created a group variable that is simply a factor variable representing each of those levels.

In my current set-up, if I want to compare differences in means between these groups, I can do something straight-forward like:

felm(outcome ~ group | team_id | 0 | 0) 

Where I cluster the standard errors at the team level. However, now I have two control conditions where I pair people up in the same way (into the four sub-groupings), but have them do different activity. I want to demonstrate that, for the same sub-group, the treatment condition has a greater mean than any of the control conditions.

In other words, I want to compare MM in the treatment group to MM in either of the two controls, and confirm that the DV of MM in the outcome group is larger than both of the control conditions.

What is the right way to do this?

My current approach is to create a condition variable that is equal to three levels, e.g. ("Main Treatment", "Control 1", "Control 2"). So treatment is a factor with four levels ("MM", "FM", "MF", and "FF"):

lm(outcome ~ treatment*condition)

But I'm not sure this is actually capturing what I intend for it to capture

---

Edits:

1. With any of the above approaches, how can I continue to do the relevant comparisons while clustering standard errors at the team level (e.g. via felm from lfe library or lm_robust from estimatr library?)

2. With the ANOVA approach described in the answer, the resulting output will look something like this:

Response: outgroup_feelings_diff
                     Df Sum Sq Mean Sq F value   Pr(>F)   
treatment             3   1276  425.23  3.8863 0.009073 **
condition             2    331  165.67  1.5141 0.220867   
treatment:condition   6   1132  188.62  1.7238 0.113025   
Residuals           592  64777  109.42

However, this doesn't tell me if a specific combination of the treatment and condition is better than some other combination -- it is presumably averaging across all combinations. Is there a way, in the ANOVA framework to ask some specific question like whether MF in the "Main Treatment" is greater than MF in "Control 1" and "Control 2." (Should I just stick with the Tukey approach if I want to do this?)

3. What is the substantive difference between the interaction approach described in #2 and creating an entirely new variable that explicitly interacts treatment and condition. For example:

treatment_condition <- paste0(treatment, condition)
lm(outcome ~ treatment_condition)

And then cycling through the relevant reference groups until I get the comparison I want.

Answer 1

Two-way anova interpretation

For the experiment you describe, you need a two-way anova (two-way since you have two factors in your experiment). The lm function (lm(outcome ~ treatment*condition)) that you describe fits this model, but there are different ways to get the results. The presentation of the results can be quite different, even though the underlying models is the same.

First way: using anova()

my_lm = lm(outcome ~ treatment*condition, data)
anova(my_lm)

The anova table will give you an estimate for each of your factors and the interaction. So this will tell you wether (1) the treatment has a significant effect on the outcome and (2) whether the condition has a significant effect on the outcome and (3) whether the treatment effect is different for certain conditions (this is the interaction effect).

Second way: using summary()

my_lm = lm(outcome ~ treatment*condition, data)
summary(my_lm)

To interpret this, we could look firstly at overall quality of the model (F-statistic and p-value, R squared and adjusted R squared. And secondly the coefficients table contains all estimated effects from the linear model. Those estimates are in the form of contrasts, meaning that there is not one estimate per 'variable', but that there is one estimate for each group (one less than number of categories, since the first category is taken as control group). On way to change the control group is using relevel().

With the same summary result, to find the estimated outcome for MM doing activity 2, you would need to sum the intercept + the estimate for MM + the estimate for activity 2.

Using this table, the estimate and the p-value for MM will tell you whether the outome of MM is on average significantly higher than the control group.

Pairwise comparisons and adding self-specified contrasts:

You want not only to compare MM to the automatically selected control group but also to other effects. One way is to get pairwise comparisons, or post hoc tests, where you compare each treatment category to each other(pairwise meaning test every possible combination for significance). To do this easily you can use the Tukey post hoc tests in TukeyHSD(). Unfortunately, they do not work with an lm object, so you need to re-fit the exact same model but with the functiom aov.

my_lm = aov(outcome ~ treatment*condition, data)
TukeyHSD(my_lm)

This will give you pairwise comparisons for each of your treatments, so basically the same as the summary output but then repeated for each group as a reference group. Also the significance values are adjusted for the fact that you do many repeated tests.

If you want to compare even more, you can take a look into contrasts. They will allow you to compare two groups of treatment categories that you specify yourself. For example you could specify a contrast for MM and FF on one side vs MF and FM on the other side.

To do this, you should look into custom contrasts. A very good explanation an be found in the answer of this question: How to set custom contrasts with lmer in R Some other reference material on contrasts here: https://rcompanion.org/rcompanion/h_01.html

And here a general reference on multi-factor anova: https://onlinecourses.science.psu.edu/stat502/node/152/

Edit:

1 - All the comparisons that I mention (even pairwise) can be specified using contrasts. I've taken a look at the felm documentation and there's also a constrast parameter that you can specify. I give you the following link that is very complete on setting up contrasts : https://rstudio-pubs-static.s3.amazonaws.com/65059_586f394d8eb84f84b1baaf56ffb6b47f.html

Basically for comparing a category A vs B and C you specify A as -1, B as 0.5 and C as 0.5 (or inverse negative / positive) as long as the sum is 0. Then it will compare the negative groups against the positive ones.

2 + 3 - For finding out whether a certain treatments works better for certain conditions, you should look at the interaction effect and see whether it is significant. Then the pairwise comparisons or specify particular constrasts for the interaction effect will detail this for particular groups.

It can become complicated here to specify those contrasts. To make it happen you could do what you propose in point 3, that is combine everything into one factor. The overall result should be the same (the F value of the omnibus test). Of course you're not going to have the split in treatment and condition and interaction effect.

Here a link that confirms this : Can you convert three-way ANOVA to one-way ANOVA?