The terms "moderation" and "interaction" are synonyms (cf., here). However, they indicate a subtle difference in how you are thinking about what you are modeling. Mathematically, neither variable included in an interaction has a different status, but often we are interested in the relationship between X and Y, and wonder if that relationship differs as a function of Z. In such cases, people sometimes refer to "moderation" instead of "interaction"; Z is the moderator here.
We usually try to visualize the relationship between two variables with a plot where we put the explanatory variable on the x-axis and the response on the y-axis. In this case, there are three variables, but you are still giving conceptual precedence to the relationship between X and Y. So it is most natural to put X on the x-axis, Y on the y-axis and then make distinct symbols or group them somehow to indicate Z. Referring to your concrete example, you are wondering about the relationship between engine displacement and fuel efficiency ("I am interested in the effect of disp
"), but wonder if that relationship varies as a function of the number of cylinders the engine has. That suggests you put disp
on the x-axis, mpg
on the y-axis, and group the symbols by cyl
.
I can demonstrate this using your example. (Note that I am taking this rather literally, which may or may not be relevant for whatever situation you really care about.) First of all, I would model the number of cylinders as a factor, not a continuous variable. Then I would use separate lines with distinct colors or types to indicate the different groupings. The existence of an interaction can be discerned by the lines not being parallel.
m = lm(mpg~cyl*disp, mtcars)
mtcars$cyl = factor(mtcars$cyl, levels=c(4,6,8), labels=c(4,6,8))
mf = lm(mpg~cyl*disp, mtcars)
AIC(m) # [1] 159.1494
AIC(mf) # [1] 153.4352
Despite the fact that using cyl
as a factor consumes more degrees of freedom, the AIC is lower for the model that does so. At any rate, here is a plot:
seq.d4 = 70:150
seq.d6 = 140:260
seq.d8 = 270:475
windows()
plot(seq.d4, predict(mf, newdata=data.frame(disp=seq.d4, cyl="4")), type="l",
xlim=c(70,475), ylim=c(5,35), xlab="disp", ylab="predicted MPG")
lines(seq.d6, predict(mf, newdata=data.frame(disp=seq.d6, cyl="6")), col="blue")
lines(seq.d8, predict(mf, newdata=data.frame(disp=seq.d8, cyl="8")), col="red")
points(mpg~disp, mtcars, subset=mtcars$cyl=="4")
points(mpg~disp, mtcars, subset=mtcars$cyl=="6", col="blue")
points(mpg~disp, mtcars, subset=mtcars$cyl=="8", col="red")
legend("topright", legend=c("4 cylinders", "6 cylinders", "8 cylinders"),
pch=1, lty=1, col=c("black", "blue", "red"))
The interpretation would be: 'Fuel efficiency generally decreases as engine displacement increases, but the relationship varies as a function of the number of cylinders. Fuel efficiency decreases more rapidly with increasing displacement when the engine has 4 cylinders than for 6 or 8 cylinders.'
If you did want to model cyl
as a continuous variable, you could use the mean of cyl
, and the mean plus/minus one SD, or at the median and quartiles of cyl
(see my answer here: How to visualize a fitted multiple regression model?).
This version adds confidence bands to the figure:
Update: It seems I misunderstood the research question the OP was trying to answer in the concrete example given.
In general, it is most natural to put the variable whose relation to the response you are interested in on the x-axis, so I have done that below. However, it is also true that the axes are better for continuous variables and grouping (e.g., by lines) is more consonant with categorical variables. Moreover, because the number of cylinders is highly correlated with the displacement, there are no 4 cylinder engines with large displacements in the dataset and the model's predicted value there is impossible. In addition, it is harder to combine the data with the predicted values coherently in this arrangement. For those reasons, the figures above might still be preferable. Nonetheless, I produce the alternate figure below. To do so, I computed the predicted values for the cylinder numbers at the mean of displacement, and plus/minus 1 SD.
md = mean(mtcars$disp); md # [1] 230.7219
mdH = mean(mtcars$disp) + sd(mtcars$disp); mdH # [1] 354.6606
mdL = mean(mtcars$disp) - sd(mtcars$disp); mdL # [1] 106.7832
ci = predict(mf, newdata=expand.grid(cyl=c("4","6","8"),disp=c(md,mdH,mdL)),
interval="confidence", level=.68)[,-1]
windows()
plot(c(4,6,8), predict(mf, newdata=data.frame(disp=md, cyl=c("4","6","8"))), type="l",
ylim=c(0,30), xlab="number of cylinders", ylab="predicted MPG", xaxp=c(4,8,2))
lines(c(4,6,8), predict(mf, newdata=data.frame(disp=mdH, cyl=c("4","6","8"))), col="blue")
lines(c(4,6,8), predict(mf, newdata=data.frame(disp=mdL, cyl=c("4","6","8"))), col="red")
arrows(x0=c(4,6,8), y0=ci[1:3,1], y1=ci[1:3,2], code=3, length=.1, angle=90)
arrows(x0=c(4,6,8), y0=ci[4:6,1], y1=ci[4:6,2], code=3, length=.1, angle=90, col="blue")
arrows(x0=c(4,6,8), y0=ci[7:9,1], y1=ci[7:9,2], code=3, length=.1, angle=90, col="red")
legend("bottomright", legend=c("mean disp", "+1 SD disp", "-1 SD disp"), lty=1,
col=c("black", "blue", "red"))