# How to test whether the relationship between an IV and a DV varies depending on levels of two other variables

Let's say I have 2 different conditions participants can be assigned to (stored in the variable condition). Within each condition, participants can be assigned to 1 of 4 treatments.

I have an outcome, Y, and a moderator variable, X and I want to know whether the relationship between X and Y depends on (and in what ways) the specific cell a participant was assigned to (treatment * condition).

My current thought about how to model this is as the following:

lm(Y ~ treatment/X -1,
data = df %>% filter(condition == "Condition A"))


And then run another model where you filter to the other condition

lm(Y ~ treatment/X -1,
data = df %>% filter(condition == "Condition B"))


This results in the following table of coefficients:

                                           Estimate Std. Error   t value   Pr(>|t|)  CI Lower  CI Upper  DF
treatment1                                3.4278052   3.337427  1.027080 0.30515602 -3.138188  9.993798 321
treatment2                               14.3758462   3.736847  3.847052 0.00014421  7.024042 21.727651 321
treatment3                                6.5085711   4.007973  1.623906 0.10537746 -1.376643 14.393785 321
treatment4                                2.3855668   3.034684  0.786101 0.43238845 -3.584815  8.355948 321
treatment1:X                             -0.0133485   0.054582 -0.244557 0.80695532 -0.120732  0.094035 321
treatment2:X                             -0.0478219   0.061340 -0.779624 0.43618685 -0.168501  0.072857 321
treatment3:X                              0.0876736   0.078806  1.112525 0.26674492 -0.067368  0.242715 321
treatment4:X                             -0.0036413   0.058255 -0.062506 0.95019843 -0.118250  0.110968 321


I assume the coefficients with the interactions measure the within-treatment effect of the independent variable on the dependent variable. Is this the appropriate way to test this kind of question? Is there a cleaner way to consolidate all the tests into a single model (instead of subsetting)?

Another thought was the following model:

lm(Y ~ treatment:condition:X -1,
data = df)


Which only models the three-way interactions between treatment, condition, and X. But I'm not exactly sure that's what I want.

I would do

lm(Y ~ X * condition * treatment, data = df)


and then look at the coefficients for the following interaction terms to see if they're significantly different from zero:

• X:condition: the average effect of X on Y depends on the type of condition
• X:treatment: the average effect of X on Y depends on the type of treatment
• X:condition:treatment: the average effect of X on Y depends on both the type of condition and the type of treatment

Note that these interpretations assume the coefficients to be significantly different from zero.

Edit after the discussion in the comments:

I've simulated some data that might resemble yours:

set.seed(123)
x <- rnorm(96, 10, 2)
treatment <- c(rep(1,24), rep(2,24), rep(3,24), rep(4,24))
condition <- rep(rep(c('A','B'), each = 12), 4)
y <- x * treatment + if_else(condition == 'A', rnorm(96, 4, 1), 0)

data <- data.frame(
x = x,
treatment = factor(treatment),
condition = factor(condition),
y = y
)

data %>%
ggplot(aes(x = x, y = y, shape = treatment, colour = condition)) +
geom_point()


Then I've fitted the linear model I mentioned above:

model <- lm(y ~ x * treatment * condition, data)


And then I've plotted the fitted regression lines on top of the data to get an idea of how the different treatments and conditions might impact the effect of X on Y:

data %>%
mutate(fitted_values = fitted(model)) %>%
ggplot(aes(x = x, y = y, shape = treatment, colour = condition)) +
geom_point() +
geom_line(aes(y = fitted_values))


• Thanks, Adria! I think the issue with this approach is that, because there are multiple treatments, the coefficient: X:treatment is actually interpreted as the average effect of X on Y, given a specific treatment, relative to the left-out base case. That is, with these fully interacted models, one treatment and one condition will be left out, so all the coefficients are just interpreted as relative to the reference case. But I want a more general answer of whether X moderates the effect of treatment on Y, for specific values of condition. Should I just rotate the base case? Sep 27 at 16:52
• @Parseltongue Could you use the fitted model to get a good idea of the effects you're interested in? I'll update my answer to include a simulated example to show what I mean. Sep 27 at 17:12
• Oh, that's a very nice idea, especially just to visualize if there is something meaningful there. Do you know if there's a more robust way to get the actual slope from all those fitted regression lines, and to test whether they're significantly different from 0? Sep 27 at 17:33
• @Parseltongue That's what the interaction terms are. For example, if the coefficient for X:condition is significantly different from zero (small p-value) then that's evidence that the impact of X on Y is different depending on the condition i.e. different slopes for the blue and red lines. Sep 27 at 17:40
• Thanks so much. Could you help me understand a perennial confusion of mine, then? Let's say the coefficient on conditionA:X is .0473 and it's not significant (p=.879) That coefficient and p-value is relative to the left out condition conditionB, correct? What it's saying is that, within condition A, there is a .0473 greater increase in the DV for every 1 unit change in X than there is in condition B. Or is it saying, just in the abstract, within condition A, a one unit increase in the X is associated with a .0473 increase in the DV? Sep 30 at 14:53