# How to Illustrate Continuous-Continuous Interactions

What is the best way to illustrate an effect modification from a general linear model? I can use GraphPad Prism or R; can anyone point me resources showcasing how to produce publication ready plots for this purpose ?

Well, I don't know about "publication ready", but you can try using the effects package in R to obtain predictor effects plots. See https://cran.r-project.org/web/packages/effects/vignettes/predictor-effects-gallery.pdf for details and also the R code below.

Example 1 [Categorical by Continuous Interaction]: Let's say you fit the model below in R and are interested in obtaining predictor effects plots from this model:

model <- lm(mpg ~ hp + cyl*wt, data = mtcars)


where cyl was converted to a factor prior to fitting the model:

mtcars$$cyl <- factor(mtcars$$cyl)


The summary of the model fit is as follows:

> summary(model)

Call:
lm(formula = mpg ~ hp + cyl * wt, data = mtcars)

Residuals:
Min     1Q Median     3Q    Max
-3.855 -1.381 -0.312  1.291  4.893

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  41.11394    3.13741  13.104 1.06e-12 ***
hp           -0.02229    0.01146  -1.945 0.063057 .
cyl6         -8.53030    8.99369  -0.948 0.351966
cyl8        -12.68744    4.85365  -2.614 0.014940 *
wt           -5.51603    1.29382  -4.263 0.000251 ***
cyl6:wt       2.27116    2.97850   0.763 0.452887
cyl8:wt       3.34995    1.54749   2.165 0.040153 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.327 on 25 degrees of freedom
Multiple R-squared:  0.8798,    Adjusted R-squared:  0.8509
F-statistic: 30.49 on 6 and 25 DF,  p-value: 2.447e-10


For this model, you can visualize the (partial) effect of hp using the code:

effect.hp <- predictorEffect("hp", model)

effect.hp

as.data.frame(effect.hp)

plot(effect.hp)


If you look at the output of as.data.frame(effect.hp):

> as.data.frame(effect.hp)
hp      fit        se    lower    upper
1  52 21.10513 1.2771874 18.47472 23.73555
2 120 19.58924 0.7274522 18.09103 21.08746
3 190 18.02876 0.8143966 16.35148 19.70604
4 260 16.46829 1.4436652 13.49500 19.44157
5 340 14.68488 2.2986623  9.95070 19.41907


you can see that - by default - the effect of hp is visualized across 5 separate hp values (ranging from 52 to 340). Behind the scenes, the effects package also chooses some "typical" values for the remaining predictor variables in the model (namely cyl and wt) in order to evaluate the desired effect of hp. You can access these values with the command:

effect.hp$model.matrix  and see that they are as follows:  (Intercept) hp cyl6 cyl8 wt cyl6:wt cyl8:wt 1 1 52 0.21875 0.4375 3.21725 0.7037734 1.407547 2 1 120 0.21875 0.4375 3.21725 0.7037734 1.407547 3 1 190 0.21875 0.4375 3.21725 0.7037734 1.407547 4 1 260 0.21875 0.4375 3.21725 0.7037734 1.407547 5 1 340 0.21875 0.4375 3.21725 0.7037734 1.407547 attr(,"assign") [1] 0 1 2 2 3 4 4 attr(,"contrasts") attr(,"contrasts")$cyl
[1] "contr.treatment"


This output makes it easy to see that R uses the observed mean value of continuous variables such as wt in its default effects calculation. Indeed:

mean(mtcars$wt)  reveals the mean of wt to be equal to 3.21725. For categorical variables, R uses the proportion of observations falling into each of the categories of that variable which was not treated as reference. The cyl variable has 3 levels - 4 cylinders, 6 cylinders and 8 cylinders - and the first of these levels was treated as reference. So we need to compute the proportion of cars with 6 or 8 cylinders in the data: prop.table(table(mtcars$wt))


which gives us this output:

>     prop.table(table(mtcars$cyl)) 4 6 8 0.34375 0.21875 0.43750  Clearly, the proportions of cars with 6 and 8 cylinders represented in the mtcars data are 0.21875 and 0.43750 and these are exactly the values reported in the cyl6 and cyl8 columns of the effect.hp$model.matrix output. Note that cyl6 and cyl8 are simply dummy variables defined as:

cyl6 = 1 if a car has 6 cylinders and 0 otherwise;

cyl8 = 1 if a car has 8 cylinders and 0 otherwise.


The effect of cyl in the model can be visualized with the R commands below in a way which makes it clear that it depends on wt (since cyl is engaged in an interaction with wt):

effect.cyl <-  predictorEffect("cyl", model)

effect.cyl

as.data.frame(effect.cyl)

plot(effect.cyl,
lines=list(multiline=TRUE),
as.table=TRUE)


The command below will produce a different visualization of the effect of cyl (which also includes measures of uncertainty):

plot(effect.cyl, lines=list(multiline=FALSE), as.table=TRUE)


The effect of wt depends on cyl and can be visualized using these R commands:

effect.wt <- predictorEffect("wt", model)

effect.wt

as.data.frame(effect.wt)

plot(effect.wt,
lines=list(multiline=TRUE),
as.table=TRUE)


To add uncertainty intervals to your visualization of the effect of wt, just use:

plot(effect.wt,
lines=list(multiline=FALSE),
as.table=TRUE,
lattice = list(layout=c(3,1)))


Example 2 [Continuous by Continuous Interaction]: Let's say you fit the model below in R and are interested in obtaining predictor effects plots from this model:

model <- lm(mpg ~ hp + disp*wt, data = mtcars)


Commands such as the ones below allow you to plot the effects of disp at pre-specified values of wt and the effects of wt at pre-specified values of disp for a typical value of hp (i.e., the mean value of hp in the data):

effect.disp <-  predictorEffect("disp", model, xlevels=list(wt = c(2.5,3,3.5)))

effect.wt <-  predictorEffect("wt", model, xlevels=list(disp = c(120,300,380)))

plot(effect.disp, as.table=TRUE, lines=list(multiline=FALSE))

plot(effect.wt, as.table=TRUE, lines=list(multiline=FALSE))


You can also show the effects of disp and wt in the same graphical window; for example:

plot(predictorEffects(model, ~ disp + wt,
xlevels=list(wt = c(2.5,3,3.5), disp = c(120,300,380))),
as.table=TRUE,
lattice = list(layout=c(1,3)))


or

plot(predictorEffects(model, ~ disp + wt,
xlevels=list(wt = c(2.5,3,3.5), disp = c(120,300,380))),
as.table=TRUE,
lattice = list(layout=c(3,1)))


The plot produced by the last command is shown below.

Typically, you calculate the conditioned glm equation for several levels of your moderator (e.g., mean - 1 sd, mean, and mean + 1 sd). This can then be plotted in a scatterplot.

I recommend using ggplot in R. To make it publication ready, you can use several themes, for instance, for APA there is a dedicated theme.

• I assume this is referred to as simple slopes analysis. Anyone have good examples of publications of this sort in social sciences research (or any field). Commented Jul 29, 2019 at 12:09
• It's a bit beyond the scope of one "simple" interaction but they are very explicit about how you should interpret and illustrate interactions: public.kenan-flagler.unc.edu/faculty/edwardsj/… Commented Jul 29, 2019 at 13:35

@Isabella' s answer is great! I wanted to add a ggplot alternative to plotting the interactions

library(effects)
library(tidyverse)

model <- lm(mpg ~ hp + disp*wt, data = mtcars)

effect.disp <-  predictorEffect("disp", model, xlevels=list(wt = c(2.5,3,3.5)))

# When turning the effect.disp object into a dataframe, we see
# that it has all the elements we want
# The "fit" is the predicted mpg
# The "lower" and "upper" are the uncertainty values we need for the ribbon

effect.disp %>% as.data.frame() %>% names()
#> [1] "disp"  "wt"    "fit"   "se"    "lower" "upper"

effect.disp %>%
as.data.frame() %>%
ggplot(aes(x = disp, y = fit))+
geom_line()+
geom_ribbon(aes(ymin = lower, ymax = upper), fill = "grey30", alpha = 0.2)+
facet_wrap(~wt)+
labs(y = "mpg")


# or if we want all lines in one plot
# [it's best if we turn the "wt" variable into a factor]

effect.disp %>%
as_tibble() %>%
ggplot(aes(x = disp, y = fit, group = factor(wt)))+
geom_line(aes(colour = factor(wt)))+
geom_ribbon(aes(ymin = lower, ymax = upper, fill = factor(wt)),
alpha = 0.2)+
labs(y = "mpg")


Created on 2019-07-31 by the reprex package (v0.3.0)

• This is so cool, Lefkios! Now we have a receipe for easier customization of effects plots. 👌 Commented Jul 31, 2019 at 10:45